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THE  USE  OF  LAGRANGE  MULTIPLIERS  FOR  ELASTIC  STRUCTURE 
WITH  ARBITRARY  BOUNDARY  CONDITIONS 
SUBJECTED  TO  THERMAL  LOADING 

By 

Un-Chin  Chai 

August,  1996 

Chairman:  Dr.  Ali  A.  Seireg 

Major  Department:     Mechanical  Engineering 

The  objective  of  this  study  is  to  develop  a  procedure 

for    generating    "closed- form"    finite    series    solutions  for 

thermal  deformation  and  stress  resulting  from  any  heat  input 

or    temperature    distribution    throughout    the    thickness  of 

multilayered    structures    of    arbitrary    shapes    and  boundary 

conditions . 

The  Rayleigh-Ritz  method  is  used  in  this  study  to 
evaluate  the  deflection  in  structures  with  boundary 
conditions  which  can  be  readily  satisfied  by  a  preselected 
deflection  function.  For  structures  with  complex  boundary 
conditions  and  shapes  including  holes  with  arbitrary 
geometry,  the  analysis  is  accomplished  by  incorporating  the 
Lagrange  Multipliers  approach  in  the  Rayleigh-Ritz  method  in 


V 


order  to  satisfy  the  constraints  at  the  boundaries  in  a 
discrete  manner. 

A  simplified  finite  difference  method  is  also  utilized 
for  rapid  evaluation  of  the  transient  temperature 
distribution  in  homogeneous  and  composite  structures  due  to 
arbitrary  heat  flux  conditions. 

In  all  the  cases  considered  for  illustration,  the  CPU 
time  required  to  reach  a  converging  solution  is  found  to  be 
orders  of  magnitude  smaller  than  that  required  to  solve  this 
class  of  problems  by  finite  element  analysis.  An  important 
advantage  of  the  proposed  approach  is  that  the  solution  is  a 
continuous  finite  series  rather  than  discrete  numerical  data 
at  the  nodes  or  elements  of  the  selected  mesh. 
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CHAPTER  1 
INTRODUCTION 

The  Thermal  Loading 

Thermal  stress  analysis  of  structures  has  been  studied 
by  investigators  for  years.  The  thermal  stress  evaluation  in 
composite  laminates  is  also  extensively  discussed  in  the 
literature.  Examples  include  simply  supported  sandwich 
panels[l],  laminated  shells [2 , 3  ,  4 ,  5]  and  fiber-reinforced 
composites [6] .  Other  examples  include  the  structural 
elements  in  nuclear  power  plants [7]  where  the  thermal- 
structural  behavior  and  performance  is  of  critical 
importance.  There  are  many  techniques  which  evaluate  the 
thermal  stress  in  various  structures  such  as  the  case  of 
beams  during  hot  dip  galvanizing [8 ]  ,  plates  with  different 
material  properties [2 , 3  ,  9 , 10 , 12  ]  and  three-dimensional  solid 
structures [3 , 11]  .  Most  of  the  previous  techniques  are  based 
on  discrete  numerical  methods  such  as  finite  element 
method (F . E .M. )  or  finite  difference  method  for  calculating 
the  temperature  and  stress  distribution.  Since  the  highest 
thermal  stress  and  deformation  in  structures  often  do  not 
take  place  at  steady  state  condition,  the  transient 
temperature  distribution [12 , 13 , 14 , 15]  at  which  the  maximum 
conditions  occur  plays  an  important  role  in  thermal  stress 
analysis . 
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In  this  study,  the  transient  thermal  analysis  is 
performed  by  a  simplified  finite  difference  method [16] . 
This  technique  combines  the  principle  of  conservation  of 
energy  and  the  theory  of  heat  penetration  distance  to 
calculate  the  temperature  field  at  any  instant  of  time.  The 
finite  difference  equations  of  incremental  control  volumes 
are  generated  by  using  the  principle  of  conservation  of 
energy.  These  control  volumes  can  be  determined  by 
utilizing  the  heat  penetration  distance.  Therefore,  two- 
dimensional  and  three-dimensional  transient  thermal  analyses 
can  always  be  possible  to  reduced  to  one-dimensional 
problems . 

Finite  Element  Analysis 

Recent  developments  in  modern  digital  computers  have 
resulted  in  extensive  activity  in  the  field  of  structural 
design  and  thermal  analysis.  Of  particular  importance  is 
the  use  of  the  finite  element  method  in  this  class  of 
problems.  It  is  a  very  flexible  and  general  method  which 
can  be  applied  to  most  of  the  elastic-plastic  and  thermal- 
elastic  problems,  even  in  the  nonlinear  domain [17, 18] .  It 
also  can  handle  complex  shapes,  loading  and  material 
properties [19] . 
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Modified  Rayleigh-Ritz  Method  by  Utilizing 
Lagrange  Multiplier 

There  are  numerous  classical  methods  which  have  been 
developed  on  the  basis  of  the  minimum  energy  principle. 
Most  of  these  methods  are  based  on  the  "principle  of  minimum 
potential  energy"  such  as  the  Rayleigh-Ritz  method.  The 
application  of  the  Rayleigh-Ritz  method  in  simple  structural 
analysis  can  be  found  in  many  references [20 , 21 , 22 ] .  These 
applications  are  generally  limited  to  simple  problems,  such 
as  a  simply- supported  beam,  or  simply- supported  plate.  Since 
the  Rayleigh-Ritz  method  requires  that  the  assumed 
deflection  function  has  to  satisfy  the  boundary  conditions, 
it  may  not  be  always  possible  to  find  a  solution  which  meets 
all  the  boundary  condition  requirements  and  minimizes  the 
assumed  function  at  the  same  time. 

In  order  to  solve  problems  with  complex  shapes  and 
arbitrary  boundary  conditions,  the  Rayleigh-Ritz  method  is 
modified  to  incorporate  the  Lagrangian  approach  for  dealing 
with  equality  constraints  [23,24,25,26].  In  this  approach, 
the  boundary  conditions  are  forced  to  be  satisfied  with  the 
use  of  a  Lagrangian  function  which  includes  the  necessary 
equality  constraints  with  undetermined  multipliers.  By 
applying  this  method,  a  finite  series  solution  can  be 
obtained  for  many  difficult  elasticity  problems  with 
arbitrary  shapes  and  boundary  conditions. 
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The  main  objective  of  this  study  is  to  developed  a 
closed-form  solution  for  thermo-elastic  problems  with 
arbitrary  shapes  and  boundary  conditions.  The  proposed 
method  is  divided  into  two  parts:  the  first  part  generates 
the  solution  for  the  transient  temperature  distribution  and 
transforms  it  to  an  equivalent  mechanical  loading,  the  other 
part  evaluates  the  deformation  and  thermal  stress  for  the 
structure  subjected  the  equivalent  mechanical  loading.  By 
combining  the  two  techniques,  the  thermo-elastic  problems 
can  be  solved  in  a  very  small  fraction  of  the  CPU  time 
required  when  using  the  F.E.A.  Furthermore,  the  solution  is 
obtained  in  the  form  of  a  finite  series  which  is  applicable 
to  the  entire  geometry  without  the  need  for  preprocessing, 
postprocessing  or  making  decisions  on  the  structure  of  the 
mesh  for  each  problem. 

Several  numerical  examples  will  be  considered  to 
evaluate  the  merits  of  the  proposed  approach.  A  detailed 
illustration  of  the  proposed  process  with  the  calculated 
numerical  data  will  be  given  in  the  Appendix.  It  can  be 
used  as  a  demonstration  of  the  different  steps  involved  in 
the  analysis  and  the  numerical  results  from  each  step. 

All  the  considered  examples  were  also  analyzed  by  using 
the  finite  elements  method  for  comparison.  The  finite 
element  models  are  created  and  analyzed  by  commercial  F.E.M. 
software  (COSMOS/M  version  1.71).  Both  the  F.E.M.  and  the 
proposed  approach  are  implemented  on  a  486-DX2  personal 
computer . 


CHAPTER  2 
THE  TEMPERATURE  DISTRIBUTION 

The  General  Formulation  of  the  Finite  Difference  Algorithm 

Before  solving  for  the  deformation  of  the  structure 
under  thermal  loading,  it  is  necessary  to  perform  the 
thermal  analysis  by  which  the  temperature  field  in  this 
structure  is  achieved.  There  are  many  methods  that  have 
been  developed  to  resolve  the  heat  transfer  problems.  In 
this  dissertation,  the  transient  heat  transfer  problems  are 
solved  by  using  the  finite  difference  method [16] .  The 
general  concepts  of  finite  difference  are  based  on  the 
principle  of  conservation  of  energy (the  first  law  of 
thermodynamics)  and  the  theory  of  heat  penetration  distance. 

From  the  principle  of  conservation  of  energy,  the 
temperature  at  each  control  volume  is  determined  by  sets  of 
difference  equations.  In  this  approach,  the  control  volumes 
have  to  be  defined  beforehand.  The  theory  of  control  volume 
is  discussed  in  later  sections  of  this  chapter  for  different 
types  of  shapes . 

After  formulating  the  difference  equations  for  the 
temperatures  of  the  control  volumes,  the  heat  penetration 
distance  at  that  instant  of  time  has  to  be  monitored.  If 
the  control  voliome  is  beyond  the  heat  penetration  distance. 


5 


6 


there  is  no  heat  transfer  on  the  boundary  of  the  control 
volume.  Therefore,  by  using  the  temperature  difference 
equations  and  the  theory  of  heat  penetrate  distance,  the 
temperature  of  each  control  volume  can  be  evaluated.  The 
temperature  field  of  the  structure  can  be  numerically 
determined  at  any  instant  by  the  calculated  temperatures  of 
control  volumes . 

In  later  sections,  the  finite  difference  results  will 
be  compared  with  the  results  of  finite  element  analysis 
(F.E.A.)  models.  Since  the  finite  difference  method  is  a 
simplified  one-dimensional  heat  transfer  for  all  kinds  of 
geometrical  shapes,  it  is  expected  that  the  finite 
difference  method  would  be  much  faster  than  F.E.A. 

One-Dimensional  Problems 

For  one-dimensional  single-material  problems.  The  heat 
flows  only  in  one  direction.  For  problems  of  beams  with 
surface  heat  input,  this  direction  is  the  direction  of  the 
thickness  of  the  beam.  In  figure  2.1,  the  thickness  of  the 
beam  is  divided  into  a  number  of  layers  whose  thickness  is 
the  thickness  of  the  control  volume.  The  volume  and  surface 
area  of  heat  transfer  of  each  control  volume  are  defined  in 
equations   (2.1)   and  (2.2). 

The  surface  area  of  each  control  volume  is 

A  =  Ai  =  A2  =  Ai=...=  An  =  An  +  \  =  bL  (2.1) 

The  volume  of  each  control  volume 


(2 .2a) 


V2  =  Vi=...=  Vn-i  =  bLAh  (2.2b) 

where  b  :  width  of  beam 
L  :  length  of  beam 
h  :  thickness  of  beam 
Ah  :  thickness  of  layers  except  first  and  last 
layers 


layer  k 


layer  n+1 


Ah/2 


(b) 


Figure  2,1 


(a)  The  beam  configuration  with  heat  source  Q 
(b)  The  control  volume  of  the  beam 
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The  assumptions  for  this  case  are  as  given  below: 

1.  The  temperature  is  assumed  uniform  within  any  layer. 

2 .  This  analysis  considers  the  transient  temperature 
variations  along  the  depth  of  the  beam  and  neglects 
the  less  significant  variations  along  the  length. 

3.  The  beam  is  subjected  to  a  uniform  heat  rate  Qin  at 
its  external  upper  surface. 

4 .  There  is  no  heat  loss  from  the  boundary  of  the  beam 
to  the  surrounding. 

The  conservation  of  energy  principal  shown  in  equation 
(2.3)   is  essentially  an  energy  balance.  The  stored  energy  in 
the  material  is  what  makes  the  temperature  of  the  material 
rise. 

Qin    —    Qout    =    Qstored  (2.3) 

The    finite    difference    scheme    is    a    simplified  one 

dimensional    finite   difference   analysis.    The   beam   shown  in 

figure  2.1  is  divided  into  n+1  layers.   The  temperatures  at 

every  surf  ace  (upper  or  lower)  are  the  main  parameters  in 
this  case.     Their  evaluation  is  based  on  the  conservation  of 

energy  principal  for  a  control  volume (the  shade  area  in 
figure  2.2). 
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^k-l,t 


^k+l,t 


<k-i  > 


4 

Heat 


Control 
volume 


Figure  2.2    The  heat  transfer  for  the  control  volume 


KA(Tk  -  1.  t  -  Tk,  t)  ,  KAt 

Qc.  k  =   At  =   A{Tk  -  1,  t  -  Tk.  t) 

Ah  Ah 


(2.4) 


^     Tk.  t  -  Tk  +  1,  t  KMt 
Oc  k  +  1  =  fCAAt   =  ——  (Tk.  t  -  Tk  +  1.  t))  (2.5) 


Ah 


Ah 


Qstored    =    pCVk{Tk,  t  +  1  -  Tk.  t) 

where  K  :  thermal  conductivity 


(2.6) 


Applying  equation  (2.3),   the  temperature  Tk,t+i  after  a 
time  step  At  can  be  presented  in  equation  (2.7)  . 


„,„                    >       KAAt  ,  JCAAt 
pcVkiTk.  t     1  -  Tk.  t)  =  — (Tk  -  i.t  -  Tk.  t)  (Tk.  t  -       +  1,  t) 


Ah 


Ah 


„                  KAtA    ,  2KAtA 
Tk.t  *  1  =   (Tk  -  1,  t  +  Tk  +  1,  t)  +  (1  —yTk.  t 


pcVkAh 


pcVkAh 


(2.7) 


Similar  expressions  can  be  obtained  for  the  temperature 
at  the  bottom  and  top  surfaces.  The  difference  equation 
(2.8)  can  be  written  as: 
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For  the  top  surface: 


Q 


in 


To,, 


Figure  2.3  Heat  transfer  for  the  first  layer 


pcVi(To.  t  +  1  -  To,  t)  =  AQinAt  - 


KAAt 


(To,  t  -  Ti,  t) 


2AtA  .  2KAtA  2KAtA 

To,  t  +  1  =   Qin  +  (1  )To,  t  +   Ti,  t  (2.8) 

pcVi  pcViAh  pcViAh 

At  the  bottom  layer  of  the  beam  shown  in  figure  2.4,  it 

is     assumed     that     the     heat     does     not     transfer     to  the 

surrounding,    and  the  only  inputs  are  from  the  upper  layer. 

The  difference  equation  can  be  written  as: 


'^n-l.t 


T 

^n,t 


Figure  2.4  Heat  transfer  for  the  last  layer 
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pcVniTn,  t  +  1  -  Tn,  t)  = 


KAAt 


{Tn   -  l.t  - 


Ah 


Tn.  t) 


Tn.  t  +  1  - 


2KAAt 


2KAAt 


pcVnAh 


Tn 


-  1. 1  +  (1  - 


pcVnAh 


yrn.  t 


(2.9) 


The  current  time  is  automatically  updated  in  the 
program  and  it  can  determine  the  distance  that  the  heat 
penetrates  into  the  beam.  This  penetration  distance  shown 
in  equation  (2.10)  is  modified  in  the  temperature  difference 
equations  based  on  the  particular  instant  of  time  according 
to  the  following  relationship: 

tp  =  0.0891^  or 


where  tp  :  penetration  time  (sec) 
Dp  :  penetration  distance 
P     :   thermal  dif f usivity (K/pc ) 

When  the  distance  between  the  heat  source  and  heat 
exchange  surface  of  the  control  volume  is  more  than  the 
penetration  distance,  no  heat  exchange  occurs.  The  general 
procedure  is  therefore  as  described  below: 

1.  Input  the  dimensions  of  the  beam: 

beam  thickness (h) 
beam  width (b) 
beam  length (L) 

2 .  Input  the  material  constants 

Density  of  the  material (p) 
Thermal  conductivity (k) 
Specific  heat(c) 


(2.10) 
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3.  Input  heat  source  (Qin)  and  total  heating  time,  t 

4.  Set  up  time  step (At) 

5 .  Calculate  the  volume  and  surface  area  of  the  control 
volume 

6 .  Record  the  current  time  =  time+  At 

7 .  Compute  how  many  layers  have  heat  exchange . 

8.  Solve  the  temperature  difference  equation 

9.  Repeat  steps  6  to  8  until  the  current  time  is  equal 
to  the  total  heating  time. 

The  flow  chart  for  the  program  is  shown  in  figure  2.5. 

For  the  case  of  a  composite  beam  shown  in  figure  2.6, 
the  surface  area  of  the  control  volume  for  each  material  and 
the  volume  of  control  voliime  for  the  material  1  are  the  same 
as  in  the  single  material  case.  Due  to  the  theory  of  heat 
penetration  distance,  the  volume  of  control  vol\ime  for  the 
rest  of  the  materials  can  be  determined  in  equation  (2.12). 

For  the  layers  of  the  first  material: 

V,,=V„,,,,,  =  ^  (2.11a) 
2 

Vk.\  =  bLiJi\,     k=2  to  nl  (2.11b) 
For  the  rest: 

Ahi  =  J-r^Ahi 


P 


1 


Vi./  =  V«2  +  i,,=— —  (2.12a) 
Vk.i  =  bLAhi,  k=2  to  n2  (2.12b) 
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Figure  2.5  The  flow  chart  of  one  dimensional  heat  transfer 
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Figure  2.6  The  configuration  of  composite  beam 
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It  is  more  convenient  to  define  the  heat  density,  HD, 
for  each  layer  in  the  different  material. 
For  a  layer  of  material  i, 

HDk  =  piCiVj.i^Jii  (2.13) 

Where  k  =  ^ni  +  j  is  the  global  count  of  the  layer  number 

The  temperature  difference  equations,  as  modified  by 
replacing  the  pcV  part  into  the  HD,  are  given  in  equations 
(2.14) ,  (2.15)  and  (2.16) . 

For  the  first  layer: 

2AtAAhi  .  2fCiAtA,  2KiAtA 

To,  t  +  1  =   Qin  +  (1  yio.  t  +   Ti,  t 

HDi  HDi  HDi 

(2.14) 

For  the  last  layer: 

2Km&tA  „  2KnAtA^ 

Tn,  t  +  1    =    Tn  -  1,  t  +  (1  )Tn,  t  (2  .  15  ) 

HDn  HDn 

where  Ahi:  the  layer  thickness  of  the  last  material 

For  the  remaining  layers  in    material  j : 

KkAtA  ,  ,      „  2KkAtA^ 

Tk,  t  +  1  =   {Tk  -  1,  t  +       +  1,  t)  +  (1  YTk.  t  (2.16) 

HDk  HDk 

As  can  be  seen,   the  general  procedure  is  the  same  as  in 

the  single  material  case.      The  only  difference  between  the 

two    cases    is    the    algorithm    used    to    generate    the  finite 

difference  equations. 
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Q 


h=0.1 


L=0.1 

(a) 


(b)  .. 
Figure  2 . 7 

(a)  The  configuration  of  the  beam  with  thermal  loading 
(b)  The  F.E.A.  model  of  Example  2.1(66  nodes,   50  elements 

and  100  time  steps) 
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Illustrative  example 
Example  2 . 1 

In  this  case,  a  homogeneous  beam  is  subjected  to  a  heat 
flux  at  its  upper  surface.  The  configuration  of  the  beam  is 
shown  in  figure  2.7a.  The  finite  element  model  is  shown  in 
figure  2.7b. 

The  material  properties,    initial  heat  input  conditions 
and  the  dimension  of  the  beam  are  given  below: 
For  Steel: 

Density  of    material (steel) , p=7900  (Kg/m^) 
Thermal  conductivity (steel) ,  k  =  45  (W/m-°C) 
Specific  heat(steel),  c=460  (J/kg-°C) 

heat  input,   Q=10  Watt 
Heat  time,   t=60  (sec) 

the  dimension  of  the  beam: 
beam  thickness,  h=0.1  (m) 
beam  width,  b=0.01  (m) 
beam  length,  L=0.1  (m) 

The   performance   of    the    finite   dif ference (F.D. )    method  is 

compared    with    the    finite    element    analysis (F.E. A. )  with 

regard  to  temperature  distribution  as  shown  in  table  2.1  and 

figure   2.8,    as   well   as    the   CPU   running   time   as    shown  in 

table    2.2.       It    can   be    seen    that   very   close    results  are 

obtained  by  the  F.D.  method  and  the  F.E. A.     The  F.D.  method, 

however,   is  much  faster  than  the  F.E. A. 


18 


Table  2 . 1  The  comparison  of  temperature  distribution 
between  F.E.A.  and  F.D. 


thickness 

F.D. 

F.E.A. 

0 

216 . 4819 

203 

0.01 

90 . 1579 

90  .3 

0.02 

50.6328 

51.9 

0.03 

29.7326 

31 . 1 

0.04 

17.3734 

18 .  6 

0 . 05 

9 . 8916 

10 . 9 

0 . 06 

5 . 4076 

6.21 

0.07 

2 .7893 

3  .46 

0.08 

1.3114 

1.91 

0.09 

0.4953 

1.14 

0.1 

0 

0.91 

ooooooooo 
oc>c><6<So<6<6c> 


thickness 


Figure  2.8  The  comparison  of  temperature  distribution 
through  the  thickness  between  F.E.A.  and  F.D. 
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Table  2 . 2     the  comparison  of  running  time  of  CPU 
between  F.E.A.  and  F.D. 


F.E.A. 

F.D. 

CPU  Time  (sec) 

56 

5 

Example  2 . 2 

This  example  considers  a  composite  structure  with  two 
materials,  Steel  and  Aluminum,  subject  to  a  heat  input  at 
the  upper  surface.  The  configuration  of  the  beam  with  the 
thermal  loading  is  shown  in  figure  2.9a.  The  finite  element 
model  is  shown  in  figure  2.9b. 

The  material  properties  and  initial  conditions  for  this 
case  as  follows: 

For  Steel: 

Density  of    material  of  Steel  ,p=7900  (Kg/m^) 
Thermal  conductivity  of  Steel,  k  =  45  (W/m-°K) 
Specific  heat  of  Steel,   c=460  (J/kg-°K) 

For  Aliaminiam: 

Density  of    material  of  Aliaminum,  p=2700  (Kg/m^) 
Thermal  conductivity  of  Aluminum,  k  =  200  {W/m-°K) 
Specific  heat  of  Aliaminum,  c=900  (J/kg-°K) 
Heat  input  time,   t=120   (sec),  Qin=l.  (Watt) 
beam  width  b=0.01  (m) 

The  performance  of  the  finite  difference  method  is 
compared  with  the  finite  element  analysis  with  regard  to  the 
temperature  distribution  as  shown  in  table  2.3  and  figure 
2.9,  as  well  as  the  CPU  running  time  as  shown  in  table  2.4. 
Again,  the  results  of  F.D.  and  F.E.A  are  very  similar  and 
the  F.D.   is  much  faster  than  the  F.E.A. 
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0.1  m 

steel 

0.1  m 

r 

Aluminium 

0.2m 

(a) 


Steel 
Aliominium 


(b) 

Figure  2 . 9 

The  configuration  of  the  beam  with  thermal  loading 
)   The  F.E.A.  model  for  Example  2.2    (12  6  nodes,  100 
elements  and  100  time  steps) 
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Table  2 . 4     the  comparison  of  running  time  of  CPU 
between  F.E.A.   and  F.D. 


F.E.A. 

F.D. 

CPU  Time 

75 

5 

oooo°oooo° 


thickness 


Figure  2.10  The  comparison  of  temperature  distribution 

between  F.E.A.  and  F.D. 


22 


Table  2.3  The  comparison  of  temperature  distribution 
between  F.E.A.   and  F.D. 


thickness 

F.D. 

F.E.A. 

0 

9.2623 

9.25 

0.01 

7  .3997 

7.39 

0.02 

5.8079 

5.8 

0.03 

4.4745 

4.46 

0.04 

3 .3805 

3  .37 

0.05 

2.5011 

2.49 

0.06 

1.8077 

1.8 

0.07 

1.2699 

1.27 

0.08 

0.857 

0.86 

0.09 

0.5397 

0.552 

0.1 

0.2886 

0.316 

0.11 

0 .2726 

0.267 

0.12 

0.2566 

0.226 

0.13 

0.2274 

0.192 

0.14 

0.1887 

0.163 

0.15 

0.1521 

0.14 

0.16 

0.1255 

0.122 

0.17 

0.0988 

0.108 

0.18 

0.0634 

0.0984 

0.19 

0.0252 

0.0928 

0.2 

0 

0.0909 

The  two-Dimensional  Problem 
For  the  case  where  the  heat  input  is  applied  only  at  a 
portion  of   the  upper   surface   of   the  beam,    the   surface  of 
control  volume  is  no  longer  one  dimensional (planar) .      It  is 
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assumed  that  all  the  points  at  the  surface  of  the  same 
control  volume  have  the  same  heat  penetration  time.  Since 
the  heat  flows  radically  from  the  position  of  heat  input, 
the  shape  of  control  volume  for  the  single  material  beam, 
shown  in  figure  2.11,   is  therefore  cylindrical. 

The  surface  area  and  volume  of  control  volume  have  to 
be  consider  in  two  parts:  left-hand-side  and  right -hand- side 
of  heat  input  due  the  probability  of  non-symmetry  in  the 
shape  with  respect  to  the  heat  input.  At  each  boundary, 
there  are  four  different  types  of  control  volumes  that  need 
different  approaches  to  obtain  their  surface  areas  and 
volumes .  The  four  types  of  control  volumes  are  shown  in 
figure  2.12. 

For  type  1 (shown  in  Figure  2.12(a)): 


arc   ab  =  — 


(2.17) 


2 


area_  oab  = 


(2.18) 


4 


For  type  2 (shown  in  figure  2.12(b)): 


arc_  ajb  =  (  Qa}r 


(2.19) 


(2 .20) 


where  Sa  :   the  angle  between  line_oc  and  line  oa 


For  type  3 (shown  in  figure  2.12(c)): 


arc_  ab  =  Qbr 


(2.21) 
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Heat  flow 


Figure  2 . 11  The  control  volumes  for  single  material  beam 


sin  2eb  +  26^  2 

area_oabc  =   (2.22) 

4 

where  Sb  :  the  angle  between  line_oa  and  line_ob 

For  type  4 (shown  in  figure  2.12(d)): 

arc_  ab  =  {Qb  -  Qa)r  (2.23) 
where  9a  :  the  angle  between  line_oc  and  line_oa 
where  9b  :  the  angle  between  line_oc  and  line_ob 

,  -      sin  (29a)  +  sin(29i)  +  2(9t  -  9a)  , 
area_ocabd  =  ^ —  ^  111^^  (2.24) 


(c)  (d) 


Figure  2.12  Four  different  types  of  C.V. 
in  the  single  material  case 


The  same  procedures  is  repeated  for  the  other  side  of 
the  beam.  The  surface  area  and  volume  of  each  control 
volume  is  the  superposition  of  L.H.S.  and  R.H.S.  The  heat 
density  (HD)  of  each  layer  is  also  determined  in  the  same 
way. 

HDk  =  pcVkAh  (2.25) 


26 


Since  the  heat  penetration  time  at  any  point  on  the 
surface  of  each  control  voliome  is  the  same,  it  is  assiimed 
that  the  heat  flow  is  always  perpendicular  to  the  surface  of 
control  volume.  Therefore  this  two-dimensional-shape  heat 
transfer  problem  can  be  simplified  to  the  one  dimensional 
heat  transfer  case.  The  temperatures  of  the  control  volumes 
are  obtained  by  using  the  same  finite  difference  equations 
(2.14) ,  (2.15)  and  (2.16) . 

For  the  composite  beam,  the  surface  of  the  control 
volume  still  remains  cylindrical  in  the  first  material  part. 
However,  due  to  the  theory  of  heat  penetration  distance,  the 
surface  of  control  voliome  in  the  rest  of  materials  is  no 
longer  cylindrical.  An  illustration  of  the  profile  of  the 
control  volume  for  the  two-material  composite  beam  is  shown 
in  figure  2.13.  The  surface  of  the  control  volume  in  the 
second  material  can  be  determined  by  the  heat  penetration 
distance. 

Due  to  the  theory  of  heat  penetration  distance,  the 
ratio  between  rj  and  R-ri  is  equal  to  the  square  root  of 
ratio  between  thermal  dif fusivities  of  material  2  and  1. 


(2 .26) 


(2  .27) 


Where  pi:  thermal  diffusivity  of  material  1 
Where  P2 :  thermal  diffusivity  of  material  2 
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Figure  2 . 13  The  control  voliome  for 
two-material  composite  beam 


It  is  more  convenient  to  represent  the  profile  of  the 
surface  in  a  cylindrical  coordinate  system.  The  point  on 
the  profile  can  be  noted  as  (r,9).  By  combining  equations 
(2.28)  and  (2.29),  the  profile  equation  is  obtained  in 
equation   (2.30) . 

r2  =  hi  CSC  (6) )  (2.28) 


r  =  n  +  r2 


(2.29) 


r  =  P  -  0csc(e)  (2.30) 
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where 


Q  =  a  -  J^)hi 


p 


Technically,  if  the  end  points  of  the  profile  are 
known,  then  the  surface  area  and  volume  of  control  volume 
are  determined  by  integrating  equation  (2.30).  The  area  and 
length  integration  are  shown  in  equations  (2.31)  and  (2.32). 
There  are  still  some  details  that  should  be  discussed. 
There  are  four  different  combinations  of  end  points  for  the 
profiles  shown  in  figure  2.14  for  the  left-hand-side 
boundary.  In  each  combination,  the  end  points  of  the 
profile  will  define  the  values  of  the  volume  and  area  of 
each  control  voliame. 


area(a,  |3)  =  j  (P  -  Q  esc  (6)  fdQ 


=  p'(P  -  a)  -  2P0[ln( 


csc(P)  -  cot  (p) 


CSC  (a)  -  cot  (a) 


)  ]  -  Q'(cot  (p)  -  cot  (a)) 
(2.31) 


curve(a,  p)  =  j  (P  -  Q  csc  (0)  )dQ 

a  (2.32) 
=  P{P  -  a)  +  Q  ln(|csc  (p)  -  cot  (a)|) 

For  type  a,  the  prof  ile  (curve  from  point  a  to  point 
b)  shown  in  figure  2.14a  is  complete.  For  type  b,  shown  in 
figure  2.14b,   the  left  end  point (point  a)   of  the  profile  is 
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on  the  left  boundary  of  the  beam  and  the  lowest  point (point 
b)  of  the  profile  is  above  the  bottom  surface  of  the  beam. 


Figure  2.14  The  profiles  of  the  control  volumes  for 
difference  boundary  effect 


For  type  c,  shown  in  figure  2.14c,  the  left  end 
point  (point  a)  of  the  profile  is  on  the  common  surface  of 
the  two  materials  and  the  lowest  point  of  the  profile  is 
below  the  bottom  surface  of  the  beam.  Therefore,  the  right 
end  point  of  the  prof  ile  (point  b)  is  on  the  bottom  surface 
of  the  beam.      For  type  d,    shown  in  figure  2.14d,    the  left 


end  point (point  a)  of  the  profile  is  on  the  left  boundary  of 
the  beam.  The  right  end  point  of  the  prof ile (point  b)  is  on 
the  bottom  surface  of  the  beam.  After  calculating  the  end 
points  for  each  type,  the  area  and  volume  for  each  control 
volume  are  determined  by  utilizing  equations  (2.31)  and 
(2.32) . 

After  defining  the  area  and  volume  of  each  control 
voliame,  the  heat  density  (HD)  for  each  control  volume  is 
calculated  in  the  same  way.  By  using  the  same  finite 
difference  equations  (2.14),  (2.15)  and  (2.16),  the 
temperature  field  of  the  beam  under  thermal  loading  is 
computed.     The  procedure  is  as  given  below: 

1.  Input    the    position    of    heat    source:    The    beam  is 
divided    into    two    parts     (left-hand-side (LHS)  and 

.    right-hand-side (RHS) )    to  calculate  the  surface  area 
and  volume  of  each  control  volume. 

2.  Input    the    layer (control    volume)    thickness    in  the 
first  material 

3 .  Calculate  the  surface  area  and  volume  of  each  layer 
of  material  1  for  the  left  hand  side  of  the  beam. 

4 .  Calculate  the  surface  area  and  volume  of  each  layer 
of  material  2  for  the  left  hand  side  of  the  beam. 

5.  Repeat   (3)   and  (4)   for  right  hand  side  of  the  beam. 

6.  Superimpose  the  results   from  the  L.H.S.    and  R.H.S. 
of  the  beam. 
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7 .  Update     the     current     time     and     heat  penetration 
distance . 

8.  Update  and  solve  the  finite  difference  equations. 

9.  Repeat   (7)   to   (8)   until  heating  time  is  over. 

For    a    beam    that     is    composed    of    more    than  two 
materials,  the  procedure  is  essentially  the  same. 

Illustrated  examples 
Example  2 . 3 

In  this  example,  a  homogeneous  structure  is  subjected 
to  the  heat  input  at  a  point  on  the  top  surface  of  the  beam. 
The  configuration  of  the  beam  and  heat  input  location  is 
shown  in  figure  2.15a.  The  finite  element  model  is  shown  in 
figure  2.15b.  The  results  are  compared  with  the  finite 
element  method  in  figure  2.16  and  table  2.5  to  illustrate 
the  accuracy  and  efficiency  of  the  two  methods. 

The  material  of  the  structure:  Steel 

The  heat  input:  Qin=10  Watt  for  60  sec 

The  dimension  of  the  structure: 

L=0.2  m 
Ll=0.1  m 
h=0.1  m 
b=0.01  m 
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(a) 


(b) 


Figure  2.15   (a)  The  configuration  of  the  beam  and 
heat  input  for  example  2.3   (b)  The  F.E.A.  model  for  example 
2.3(231  nodes,   200  elements  and  100  time  steps) 


33 


250 


Thickness 


Figure  2.16  The  comparison  between  F.E.A.   and  F.D.  for 
the  temperature  distribution  through  the  thickness 
below  the  location  of  the  heat  input 


Table  2.5  The  comparison  between  F.E.A.  and  F.D. 
for  the  CPU  running  time 


F.E.A. 

F.D. 

CPU  Time 

125 

6 

Example  2 . 4 

In  this  example,   a  composite  beam  is  subjected  to  heat 
input    at    a    point    on    the    top    surface    of    the    beam.  The 
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configuration  of  the  beam  and  heat  input  location  is  shown 
in  figure  2.17a.  The  finite  element  model  with  initial 
condition  is  shown  in  figure  2.17b.  The  results  are 
compared  with  the  finite  element  method  to  illustrate  the 
accuracy  and  computational  efficiency  of  the  proposed 
method. 

The  heat  input:  Qin=10  Watt  for  120  sec 

The  dimension  of  the  structure: 

L=0.2  m 

Ll=0.1  m 

hi=0.1  m 

h2=0.1  m 

b=0.01 


LI  W  Qin 


Steel 


Aluminium 


(a) 
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Steel 
I     I  Aluminium 


(b) 

Figure  2.17(a)  The  configuration  of  the  beam  and  heat  input 
for  Ex.   2.4   (b)   The  F.E.A.  model  for  example  2.4   (441  nodes, 
400  elements  and  100  time  steps) 


300 


Thickness 


Figure  2.18  The  comparison  between  F.E.A.  and  F.D. 
for  the  temperature  distribution  through  the  thickness  the 
below  location  of  the  heat  input 
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Table  2.6  The  comparison  between  F.E.A.  and  F.D. 
for  the  CPU  running  time 


F.E.A. 

F.D. 

CPU  Time 

273 

6 

X(at  Y=0.) 


Figure  2.19  The  comparison  between  F.E.A.  and  F.D. 
for  the  temperature  on  the  top  surface 
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Three-Dimensional  Problems 
If  the  heat  input  is  on  a  portion  of  the  top  surface  of 
a  plate,  the  shape  of  control  volume  is  now  expanded  into 
three-dimensional  space.  The  surface  area  and  volume  of  each 
control  volume  of  the  structure  are  defined  by  the  principle 
of  heat  penetration  distance  as  in  the  case  of  the  one  or 
two  dimensional  problems.  Although  the  computation  becomes 
more  complex  in  three  dimensional  problems,  the  approach  is 
very  similar  to  two  dimensional  problems. 

Two  examples  are  considered  for  illustration.  One  is  a 
single-material  plate  and  the  other  is  a  composite  plate. 
Both  examples  are  subjected  to  heat  input  on  the  top 
surfaces  of  the  plate. 

Illustrate  examples 
Example  2  . 5 

In  this  example,  a  homogeneous  square  plate  is 
subjected  to  heat  input  at  the  center  of  the  top  surface  of 
the  plate.  The  plate  configuration  is  as  shown  in  figure 
2.20a.  The  finite  element  model  is  shown  in  figure  2.20b. 
Because  of  the  symmetry  of  the  geometry,  loading,  and 
boundary  conditions  of  this  problem,  only  one  quarter  of  the 
structure  is  used  for  the  F.E.A.  model.  Its  result  is 
compared  with  finite  the  element  method  in  figure  2.21  and 
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table    2.7    to    illustrate    the    accuracy    and  computational 

efficiency  of  the  proposed  method. 

The  dimension  of  the  plate: 
Li=0.2  m 
L2=0.2  ra 
Lii=L2i=0 . 1  m 
h=0.1  m 


(b) 

Figure  2.20   (a)  The  plate  configuration  and 
heat  input  Qin=10  Watt   (b)  F.E.A.  model  for  example  2.5 
(1331  nodes,   1000  elements  and  50  time  steps) 


3? 


Table  2.7  The  comparison  between  F.E.A.  and  F.D. 
for  the  CPU  running  time 


F.E.A. 

F.D. 

CPU  Time 

584 

6 

F.E.A. 
F.D. 


Thickness 


Figure  2.21  The  comparison  between  F.E.A.   and  F.D. 
for  the  temperature  distribution  through  the  thickness 
below  the  location  of  the  heat  input 


Example  2 . 6  ■  • 

In  this  case,  a  composite  square  plate  is  subjected  to 
heat  input  at  the  center  of  the  top  surface  of  the  plate. 
The  plate  configuration  is  shown  in  figure  2.22a.  The 
finite  element  model  is  shown  in  figure  2.22b.  The  result 
is  compared  with  the  finite  element  method  to  illustrate  the 
accuracy  and  computational  efficiency  of  the  proposed 
method. 
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(b) 

Figure  2.22    (a)   The  plate  configuration  and  heat  input  of 
example  2.6   (Lii=L2i=0 . 1,   Li=L2  =  0.2,   hi=h2=0.1,   Qin=10  Watt)  (b) 
F.E.A.  model  for  example  2.6(2531  nodes,   2000  elements  and 

3  0  time  steps) 
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Table  2.8  The  comparison  between  F.E.A.   and  F.D. 
for  the  CPU  running  time 


F.E.A. 

F.D. 

CPU  Time 

824 

7 

Thickness 


Figure  2.23  The  comparison  between  F.E.A.  and  F.D. 
for  the  temperature  distribution  through 
the  thickness  at  the  center  of  the  beam 


CHAPTER  3 

THERMAL  DEFLECTION  AND  STRESS  DISTRIBUTION      \     ' . 

The  General  Concept  '         '  -' 

When  a  structure  is  heated  above  the  ambient 
temperature,  it  will  expand.  Due  to  the  heat  input  into  the 
structure,  the  thermal  loading (temperature  field)  will  cause 
deflection  as  well.  The  procedure  of  obtaining  the 
temperature  field  of  the  structure  is  discussed  in  chapter 
2.  In  this  chapter,  the  main  goal  is  to  obtain  the  thermal 
deflection  and  stress  by  assuming  that  the  temperature  field 
in  the  structure  is  known. 

The  purpose  of  this  study  is  to  develop  a  new  method  to 
obtain  the  deflection  of  the  structure  under  thermal 
loading.  In  the  proposed  method,  the  thermal  loading  of  the 
structure  is  transformed  to  an  equivalent  mechanical  loading 
that  can  produce  exactly  the  same  deflection  as  the  thermal 
loading.  The  details  of  how  to  transform  the  thermal 
loading ( temperature  field)  to  an  equivalent  mechanical 
loading  will  be  discussed  later  in  this  chapter. 

After  obtaining  the  equivalent  mechanical  loading  of 
the  structure,  the  deflection  due  to  the  equivalent 
mechanical  loading  is  evaluated  by  utilizing  the  energy 
method  with  Lagrange  multipliers.      The  application   of  the 
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energy  method  with  Lagrange  multipliers  in  structural 
analysis  is  presented  in  references [23 , 24 , 25 , 26] .  The 
traditional  energy  method  , such  as  Rayleigh-Ritz  method,  is 
generally  limited  to  the  simplest  problems  where  all  the 
boundary  conditions  can  be  satisfied  by  the  assumed 
deflection     function.  Through     the     use     of  Lagrange 

multipliers,  the  problems  with  complex  boundary  conditions 
can  be  solved. 

In  this  method,  the  potential  energy  function  is 
obtained  by  assuming  a  deflection  function  composed  of  a  set 
of  coordinate  functions  On  with  unknown  coefficients.  The 
Lagrange  multipliers  are  then  utilized  to  insure  the 
satisfaction  of  any  boundary  conditions  which  are  not 
satisfied  by  the  assumed  deflection  function.  The  unknown 
coefficients  of  the  assumed  deflection  function  can  be 
determined  by  minimizing  the  Lagrangian  function  which 
includes  the  potential  energy  and  the  constraints. 

The  One-Dimensional (Beam)  Problem 
A  beam  is  first  considered  for  illustrating  the 
proposed  method.  For  example,  a  multiple  material  beam  as 
shown  in  figure  3.1  is  assumed  to  be  subjected  to  a 
temperature  distribution  through  the  thickness  as  a  result 
of  thermal  loading.  Assuming  that  the  cross  section  of  the 
beam  remains  plane,  the  neutral  axis  of  the  beam  can  be 
determined  by  geometric  consideration  as  follows: 
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Figure  3 . 1  The  cross  section  of  the  composite  beam 


The  i-layer  of  the  beam  has  a  thickness  hi.  The  yi-i 
and  Yi  are  the  distances  from  the  neutral  axis  to  the  top 
surface  and  the  lower  surface  of  i-layer.  The  yi  can  be 
expressed  as: 

i-l 

yi  =  yo  +  ^  hj  for  j=l  to  n  (3.1) 

Therefore,  for  any  layer  in  the  beam,  the  distances  of 
the  top  and  lower  surfaces  are  functions  of  the  distance 
from  the  neutral  axis  to  the  top  surface  of  the  first 
layer (yo).  Since  the  first  moment  calculated  about  the 
neutral  axis  is  zero,  the  yo  is  calculated  by  using  equation 
(3.2).  The  position  of  N.A.  can  be  therefore  determined 
from: 

n  y 

^EijydA  =  Q  (3.2) 
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The  total  flexural  rigidity  D  can  be  obtained  as  the 
sum  of  the  rigidities  of  the  different  layers  with  respect 
to  the  neutral  axis (N. A.). 

D=  2  EJ*  (3.3) 
k  =  \ 

where  I*=^+bhkdk^ 

dk=  Vk-  1  hk 

2 

The  longitudinal  thermal  expansion  of  the  beam  aT  will 
generally  cause  uniform  expansion  as  well  as  bending.  Since 
the  beam  is  free  to  expand  in  the  x-direction  and  there  is 
no  other  mechanical  loading  other  than  the  thermal  loading, 
the   resultant    stress    is    the    superposition   of    the  thermal 

stress  (CTt)  ,  uniform  stress  (Cave)  and  bending  stress  (am). 
In  order  to  satisfy  the  equilibrium  conditions,  the  total 
force  due  to  the  stress  distribution  in  x-direction  has  to 
be  zero  and  the  total  moment  about  the  neutral  axis  should 
be  also  zero. 

Therefore,   for  the  single  material  beam  with  depth  h: 

o,  =  -aET 

hll 

(5ay.=  joETdy  (3.4) 

-h/2 

am  =  ^  \aETydy 
h  -Ml 

where :       E  :  Young ' s  Modulus 

a  :  Thermal  expansion  coefficient 
T  :   Temperature  rise  value  at  any  point 
through  the  thickness  of  the  beam 
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Similarly  for  the  composite (n-layers)  beam: 


For  layer  i : 


0,(0  =  -€t/ET 


(3.5) 


n  yj 


(3.6) 


(3.7) 


The  resultant  stress   for  any  layer  of  the  beam  can  be 
calculated  as : 


The  deflection  of  the  beam  is  the  superposition  of 
uniform  expansion  and  pure  bending.  The  uniform  expansion 
is  caused  by  Oave  and  pure  bending  is  caused  by  On,.  Assuming 
that  the  length  of  the  beam  is  much  greater  than  the  beam's 
total  thickness,  the  effect  of  the  uniform  expansion  on 
deflection  can  be  ignored.  Therefore  the  deflection  of  the 
beam  with  thermal  loading  can  be  simulated  as  the  deflection 
of  the  beam  with  an  equivalent  thermal  moment  Mt  through  the 
length  of  the  beam.  This  equivalent  moment,  shown  in 
Equation  (3.9),  is  the  total  moment  caused  by  the  bending 
stress  (Gn,)  about  the  neutral  axis.  The  deflection  of  the 
beam  with  T(x,y)  can  be  simulated  as  loading  with  an 
equivalent  moment  Mt(x)   acting  on  the  beam. 


(3.8) 
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Mt(x)  =  X  ^  J  (^^^^T(x,  y)ydy  (3.9) 

J  =  l  yj-i 

Where    b:   the  width  of  the  beam 

If  the  temperature  distribution  T  is  dependent  on  y 
only,  the  Mt  is  constant  along  the  length  of  the  beam.  If 
the  T  is  a  function  of  x(the  length  of  the  beam)  and  y(the 
thickness  of  the  beam) ,  the  equivalent  thermal  moment  Mc  is 
function  of  x. 

For  the  beam  case,  the  potential  energy  is  the  sxjan  of 
strain  energy  and  the  work  done  by  loading.  The  strain 
energy  can  be  expressed  as  function  of  deflection (v)  as 
given  in  equation  (3.10).  The  work  done  by  external  force 
or  moment  also  can  be  presented  as  a  function  of  deflection. 
Therefore,  the  potential  energy  is  a  function  of  deflection. 
The  deflection  of  the  beam  v(x)  is  assumed  to  be  a  linear 
combination  of  a  set  of  prescribed  coordinate  functions 
4>n(x)   as  in  equation  (3.11). 

N 

V   =    ^  AnOnfx) 

n=l 

where  An  are  unknown  parameters 

From  the  previous  discussion,  the  thermal  loading 
T(x,y)    can    be    transformed    to    an    equivalent    moment  Mt(x). 


(3.10) 
(3.11) 
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The  work  done  by  the  thermal  loading  can  be  expressed  as 
follows : 

1 

w  =  j  Mt{v")dx  (3  .12) 

0 

If  the  beam  is  subject  to  the  constraints: 
Gj{Ai,  A2,  ■■•  ,  An)  =  0       j=l,2,...,Ic  (3.13) 
By  using  the  constraints  with  Lagrange  Multipliers,  the 
Lagrangian  function (L)   can  be  written  as  follows: 

Ic 

L  =  U  -  w  +       '^^^^  (3  . 14) 

i  =  l 

By  minimizing  the  Lagrangian  function  as  shown  in 
equations  (3.15)  and  (3.16),  the  unknown  parameters (An, 
n=l...N)  can  be  readily  evaluated.  The  deflection  of  the 
beam  is  therefore  determined  by  substituting  these 
parameters  in  equation  (3.11). 

3l 

—  =0  for  n=l,2, . . . ,N  (3.15) 
aAn 

.  .  . 

—  =0  for  3=1,2,..,   Ic  (3.16) 

The  procedure  of  the  proposed  method  is   summarized  as 
below: 

1.  Input     material     properties     and     solve      for  the 
position  of  N.A.  and  the  flexural  rigidity  D 

2.  Input  the  thermal  loading (temperature  field) 

3.  Transfer    the    thermal    loading    to     the  equivalent 
mechanical  loading 
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4.  Assume  the  series  of  deflection  function 

5 .  Calculate  the  potential  energy 

6.  Add  any  unsatisfied  boundary  conditions  with 
Lagrange  multipliers  to  the  equation  of  potential 
energy  to  form  the  Lagrangian  function 

7.  Develop  and  solve  the  set  of  linear  equations 
obtained  by  partial  derivatives  of  the  Lagrangian 
function  with  respect  to  the  unknown  coefficients 
and  the  Lagrange  Multipliers. 

8 .  Substitute  back  the  unknown  coefficients  and  the 
deflection  is  determined 

Illustrated  examples 

In  order  to  illustrate  the  accuracy  and  the 
computational  efficiency  of  the  proposed  method,  the  results 
from  the  proposed  method  are  compared  with  the  results  from 
the  finite  element  method  for  all  the  considered  examples. 

Examples  3-1  to  3-3  are  for  single  material  beams 
subjected  to  three  different  sets  of  boundary  conditions. 
Examples  3-4  to  3-6  are  two-material  composite  beams 
subjected  to  the  same  sets  of  boundary  conditions  as 
examples  3-1  to  3-3.  Examples  3-7  to  3-9  are  two-material 
composite  beams  symmetrically  structured  and  also  subjected 
to  the  same  sets  of  boundary  conditions  as  examples  3-1  to 
3-3  . 
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The  temperature  loading  for  these  nine  cases  are 
assumed  to  be  the  same.  The  temperature  difference  between 
top  and  bottom  surfaces  is  120°  and  the  temperature 
distribution  through  the  thickness  of  the  beam  is  assumed  to 
be  linear.  In  all  examples,  the  length  and  deflection  units 
are  in  meters  and  the  temperature  units  are  in  °C. 

Example  3 . 1 


L=20 


—  b=l  — 


Figure  3.2  The  configuration  of  the  beam  subjected  to 
thermal  loading  in  example  3 . 1 


For  the  homogeneous  beam  shown  in  figure  3.2,  the 
boundary  conditions  are  zero  deflection  at  both  ends  and  at 
the  middle  of  the  beam.  Since  sine  functions  already 
satisfy  these  boundary  conditions  at  both  ends,  only  one 
constraint  for  the  middle  point  is  necessary  for  this  case. 
The  assumed  deflection  function (sine  series)  is  given  in 
equation  (3.17).  The  potential  energy  of  the  beam  can  be 
expressed  by  equation  (3.18)  as  derived  from  the  assumed 
deflection  function.  The  boundary  conditions  of  the  beam  can 
be  considered  as  constraints  at  discrete  locations  as  shown 
in  equation  (3.19).  The  equilibriiim  of  the  beam  is  satisfied 
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when  the  beam  has  minimum  potential  energy.  By  using  the 
method  of  Lagrange  multipliers,  the  Lagrangian  function  is 
represented  as  equation  (3.20).  The  partial  derivatives  with 
respect    to  and   X,j    equal    zero    become    a    set    of  linear 

equations  (  (3  .  21)   and   (3.22)).     The  unknown  coef  f  icients  (A-^) 
of     the     finite     series     are     obtained     by     solving  these 
equations.      The  process   can  be   repeated  with  expansion  of 
the   series   until    the   solution   converges   with  satisfactory 
accuracy. 


An  sin  (  ))  (3.17) 

n  =  l  ^ 

Since      sin(  )sin(  )  =  0        for  m^n 

J"        L  L 

Therefore 

\\v"fdx  =  \\f  innf  An  sin  {—)){f{nKf  An  sin  (—)) 

•10  Jo  J ,  T 

n  =  l  n  =  l  ^ 

fV?dx  =  f  fn%^A.^sin^(i^)) 

Jo  Jo  T, 

n  =  l  J-i 

The  potential  energy  n  in  this  beam: 

Ue  -  W  =  -  [\v"fdK  -  \  Mt{v")dx 
2  Jo  J 

N       ^  N  (3.18) 

2     i.— n\^AJ)  -  2  ^MtiniAn 

J3  =  l,2,3...     ^  n  =  l,3,5... 

The  constraint: 

G  =  2,  An  sin(—)  =  0  (3.19) 

n  =  l  J-i 


52 


The  Lagrangian  function  is  : 
L  =  U  -  w  +  XG 


f  =  0 

dAi 

3l 


for  i=l  to  n 


(3.20) 
(3.21) 


=  0 


(3 .22) 


The  temperature  distribution  T(y)   is  as  follow: 


T{y)  =  120(1  -  ^  ^  ) 


(3.23) 


The  model  for  the  finite  element  analysis  used  for 
comparison  is  shown  in  figure  3.3.  Since  the  boundary 
conditions  and  the  thermal  loading  are  symmetrical  about  the 
center  location,  this  model  only  contains  half  of  the  beam 
with  55  nodes  and  40  elements. 


> 
> 


10 


Figure  3.3  The  F.E.A.  model  for  example  3.1 


The  results  of  the  proposal  method  for  different  number 
of  terms  N  are  shown  in  figure  3.4.  The  assumed  deflection 
function  converge  at  N=7 .  The  comparison  between  the 
proposed  method  and  the  F.E.A.  is  shown  in  figure  3.5.  The 
difference  in  the  maximum  deflection (x=16 . )  is  about  9%. 
The  CPU  time  for  calculating  deflection  by  the  F.E.A.    is  9 
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sec.  The  proposed  method  gives  the  solution  in  less  than  1 
sec . 

The  results  illustrate  that  the  proposed  method  is  much 
faster  than  the  F.E.A.  and  the  two  methods  give  very  close 
deflection  values. 


X 


5  10  15  20 


Figure  3.4  the  computed  deflections  for  N=3,5,7,9 


X 


5  10  15  20 

Figure  3.5  Deflections  by  the  proposed  method (N=7) 

and  the  F.E.A. 
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Example  3 . 2 

Consider  a  clamped- supported  homogeneous  beam  as  shown 
in  figure  3.6  subjected  to  the  same  temperature  distribution 
T(x,y)   given  in  equation  (3.23). 


1 


L=20 


b=l 

( 

h=2J 

Figure  3 . 6  The  configuration  of  the  beam  in  example  3 . 2 


By  choosing  the  same  deflection  function,  the 
formulation  of  the  total  potential  energy  is  the  same  as 
that  in  example  3-1.  The  boundary  conditions,  however,  are 
not  all  satisfied  by  this  function.  The  unsatisfied 
boundary  condition  is  v'=0  at  x=0,  which  can  be  written  as: 

G  =  ^Ann  =  0  (3.24) 

n  =  l 

The  same  procedure  as  in  example  3-1  is  applied  to 
solve  for  the  A„  and  X.  The  model  of  the  finite  element 
analysis  used  for  comparison  is  shown  in  Figure  3.7.  Since 
the  boundary  conditions  and  the  thermal  loading  are  not 
symmetrical  in  x  direction,  this  model  contains  the  full 
configuration  of  the  beam  with  105  nodes  and  80  elements. 
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Figure  3.7  The  F.E.A.  model  for  example  3.2 


V 


Figure  3.8  the  computed  deflections  for  N=5,10,20 

The  results  from  the  proposed  method  for  different 
terms  N  are  shown  in  Figure  3.8.  The  assumed  deflection 
function  converges  at  N=10.  The  comparison  between  the 
proposed  method  and  the  F.E.A.  is  shown  in  figure  3.9.  The 
difference  in  the  maximum  deflection (x=14)  is  approximately 
0.88%.  The  CPU  time  for  calculating  deflection  by  the 
F.E.A.  is  11  sec.  The  proposed  method  requires  less  than  1 
sec . 
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Again,  the  results  illustrate  that  the  proposed  method 
is  much  faster  than  the  F.E.A.  and  that  the  two  methods  give 
almost  identical  results. 


5  10  15  20 

Figure  3.9  Deflections  by  the  proposed  method (N=10) 

and  the  F.E.A. 

Example  3 . 3 

Under  certain  conditions,  for  example  a  homogeneous 
cantilever  beam,  it  is  necessary  to  use  a  more  complete 
Fourier  series,  such  as  a  series  including  both  sine  and 
cosine  functions.  The  nature  of  the  sine  functions  alone 
makes  it  difficult  to  provide  a  deflection  at  the  free  end. 

Assuming  that  the  deflection  is  a  complete  Fourier 
series  as  follows: 
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V  =  Ao  +  2^  An  COS  (  )  +  2^  Bn  sin  (  ) )  (3.25) 

The  strain  energy  stored  in  the  beam  for  this  assumed 
deflection  is: 


Ln  =  l 


w    w              4n^n3^  ^  ^  -26) 

y  y  iAnBn  —  ; —  ( COS  (rOT )  COS  (mn; )-!)}) 

m*n 

The  work  done  by  the  equivalent  moment  Mt(x)  is  still 
the  same  as  equation  (3.12).  The  boundary  conditions  are 
y=0  and  y'=0  at  x=0,  and  y"  =  0  and  y"'=0  at  x=L,  none  of 
which  are  satisfied  by  equation  (3.25).  The  constraints  are 
thus : 

N 

Gi  =  Ao  +  ^  An  =  0  (3.27) 

n  =  l 

G2  =  ^nBn  =  0  (3.28) 

n  =  l 
W 

G3  =  ^  AnH^  COS  (riTT)  =  0  (3.29) 

n  =  l 
N 

Gi  =  ^Bnn^  cos  (nn)  =  0  (3.30) 
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The  Lagrangian  function  L  can  be  formulated  as: 

4 

L  =  U  -  w  +        XiGi  (3.31) 

i  =  l 

Equating  the  partial  derivatives  of  the  Lagrangian 
function  with  respect  to  parameters  "ki  (1=1,2,3,4)  Ao,  An  and 
Bn  (n=l,2,...N)  to  zero,  a  set  of  (2N+5)  linear  equations 
with  (2N+5)  unknowns  is  obtained.  Therefore  all  the  unknown 
coefficients  can  be  solved  and  the  deflection  function  of 
the  beam  is  determined. 

The  F.E.A.  model  used  for  comparison  is  the  same  as  in 
example  3-1-2  except  for  the  boundary  conditions.  The 
results  from  the  proposed  method  for  different  number  of 
coefficients  N  are  shown  in  figure  3.10.  The  assumed 
deflection  function  converges  at  N=8 .  The  comparison 
between  the  proposed  method  and  the  F.E.A.  is  shown  in 
Figure  3.11.  The  difference  in  the  maximum  deflection (x=20 ) 
is  approximately  1.48%.  The  CPU  time  for  calculating 
deflection  by  the  F.E.A.  is  12  sec.  The  proposed  method 
only  requires  less  than  2  sec . 

Again,  the  results  illustrate  that  the  proposed  method 
is  much  faster  than  the  F.E.A.  and  that  the  two  methods  give 
almost  identical  deflections  values. 
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Figure  3.10  the  computed  deflections  for  N=6,8,10 
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Figure  3.11  Deflections  by  the  proposed  method (N=8) 

and  the  F.E.A. 
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Example  3 . 4 


L=30 


A" 


A 


1 

b=l 

=2 

1 

J 

h2= 

1^ 

Figure  3 . 12  The  configuration  and  boundary  conditions  of  the 

beam  in  example  3-4 


I. A. 

yo 

Figure  3.13  The  cross  section  of  the  beam  in  example  3-4 


For    a    composite     beam     shown     in     figure     3.12,  the 
temperature  distribution  of  the  beam  is  assumed  to  be 

T(y)  =  120(1  -  ^  t        )  (3.32) 


n  y 

J^Ei  j  ydA  = 


0 


i  =  l 


yn-  1 


yo  +  hi  yo  +  hi  +  h2  (3.33) 

j  EiydA  +      J  E2ydA  =  0 

yo  yo  +  hi 

where  hi  :  the  thickness  of  material  1( Steel) 
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h2 
El 
E2 


the  thickness  of  material  2 (Aluminum) 
the  Modulus  of  material  1 
the  Modulus  of  material  2 


The  position  of  the  N.A.  (Yo)  shown  in  figure  3.13  is 
calculated  by  using  equation  (3.33).  The  flexural  rigidity 
D     is  calculated  as  follows: 

D  =  t{Ei(—  +  hi(—  +  yo)']  +  i3[B2(        +  h2(hi  +  —  +  yof]  ] 

12  2  12  2 

where  b:  width  of  the  beam  (3.34) 

The  equivalent  moment  Mt(x)   is  obtained  as  follow: 

Mt(x)  =    j  BiaiT(x,  y)ydy  +    j  E2a2T(x,  y)ydy  (3.35) 

yo  yo  +  hi 

The  work  done  the  equivalent  moment  Mt(x)  is: 

L 

w  =  j  Mtv"dx  (3.36) 

0 

The  potential  energy  n  in  this  beam  and  the  constraint 
are  the  same  as  equation  (3.18)  and  (3.19).  The  Lagrangian 
function  L  is  also  the  same  as  equation  (3.20).  By  the  same 
procedure  as  equation  (3.21)  and  (3.22),  the  deflection 
function  of  the  beam  is  determined. 

The  F.E.A.  model  used  for  comparison  is  shown  in  figure 
3.14.  Because  of  the  symmetrical  geometry  and  boundary 
conditions,  this  model  only  considers  half  of  the  beam  that 
contains  77  nodes  and  60  elements.  The  results  from  the 
proposed  method  for  different  number  of  terins  N  are  shown  in 
Figure   3.15.      The  assumed  deflection   function  converges  at 
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N=7 .  The  comparison  between  the  results  from  the  proposed 
method  and  the  F.E.A.  is  shown  in  Figure  3.16.  The 
difference  in  the  maximum  deflection (x=26)  is  approximately 
3.10%.  The  stress  distributions  through  the  thickness  from 
the  F.E.A.  and  the  proposed  method  are  shown  in  figure  3.17. 
The  shear  stress  at  the  interface  surface  of  the  two 
materials  as  obtained  by  the  proposed  method  is  37.28  MPa. 
The  CPU  time  for  calculating  the  deflection  by  the  F.E.A.  is 
12  sec.     The  proposed  method  requires  less  than  2  sec. 

The  results  illustrate  that  the  proposed  method  is  much 
faster  than  the  F.E.A.  and  that  the  two  methods  give  very 
close  deflection  values. 


Figure  3.14  The  F.E.A.  model  for  example  3.4 


63 


Figure  3 . 16 
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Example  3 . 5 

Consider  a  clamped- supported  composite  beam  as  shown  in 
Figure  3.18  subjected  to  the  temperature  distribution  T(x,y) 
as    shown   in  equation    (3.23).      This   composite  beam   is  the 
same  as  the  beam  in  example  3.4. 

   L=30   - 


Figure  3.18  The  configuration  of  the  beam  in  example  3.5 
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Following  the  same  procedure  as  in  example  3.4,  the  Mt 
can  be  obtained.  By  choosing  the  same  deflection  function, 
the  formulation  of  the  total  potential  energy  is  the  same  as 
that  in  example  3.2.  Since  the  boundary  conditions  are  not 
all  satisfied,  the  unsatisfied  boundary  condition  (v'=0  at 
x=0)   can  be  written  as: 


G  =  ]^  Ann  =  0 


n  =  l 


The  same  procedure  as  in  example  3.2  is  applied  to 
solve  for  the  An  and  X.  The  F.E.A.  model  used  for  comparison 
is  shown  in  Figure  3.19.  Since  the  boundary  conditions  and 
the  thermal  loading  are  not  symmetrical  in  the  x  direction. 
This  model  contains  the  full  configuration  of  the  beam  with 
147  nodes  and  120  elements. 


^  A 
Figure  3.19  The  F.E.A.  model  for  example  3.5 
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The  results  from  the  proposed  method  for  different 
number  of  term  N  are  shown  in  figure  3.20.  The  assumed 
deflection  function  converges  at  N=10.  The  comparison 
between  the  proposed  method  and  the  F.E.A.  is  shown  in 
Figure  3.21.  The  difference  in  the  maximum  deflection (x=20 ) 
is  approximately  1.69%.  Since  the  temperature  loading  is 
the  same  as  in  example  3.4,  the  stress  distribution  through 
the  thickness  of  the  beam  and  slip  stress  is  the  same  as 
that  in  example  3-4.  The  CPU  time  for  calculating 
deflection  by  the  F.E.A.  is  11  sec.  The  proposal  method 
requires  less  than  1  sec. 
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Again,  the  results  illustrate  that  the  proposed  method 
is  much  faster  than  the  F.E.A.  and  that  the  two  methods  give 
almost  identical  deflection  values. 

Example  3 . 6 

In  this  example,  we  consider  a  composite  beam  with  the 
same  boundary  condition  as  in  example  3.3  and  the  same 
configuration  of  the  beam  as  in  example  3.4.  The  beam  in 
this  case  is  shown  in  figure  3.22. 


68 


  L=30  

Figure  3.22  The  configuration  of  the  beam  in  example  3-6 


Following  the  same  procedure  as  in  examples  3.3  and 
3.4,  the  deflection  is  obtained.  The  F.E.A.  model  used  is 
the  same  as  in  example  3.5  with  different  boundary- 
conditions.  The  results  from  the  proposal  method  for 
different  number  of  terms  N  are  shown  in  figure  3.23.  The 
assumed  deflection  function  converges  at  N=8 .  The 
comparison  between  the  proposed  method  and  the  F.E.A.  is 
shown  in  Figure  3.24.  The  difference  in  the  maximum 
def lection (x=30)  is  approximately  3.09%.  The  CPU  time  for 
calculating  deflection  by  the  F.E.A.  is  12  sec.  The 
proposed  method  requires  less  than  2  sec. 

Again,  the  results  illustrate  that  the  proposed  method 
is  much  faster  than  the  F.E.A.  and  those  two  methods  give 
almost  identical  deflection  values 


Figure  3.24  Deflections  by  the  proposed  method (N=8) 

and  the  F.E.A. 
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Example  3 . 7 


L=30 


hi=l 


h2=l 


£ 


T 


h3  =  l-' 


Figure  3.25  the  configuration  and  boundary  conditions  of  the 

beam  in  example  3 . 7 


For    the    composite    beam    shown    in    figure    3.25    and  a 
temperature  distribution  as  given  by  equation  (3.37): 


T(y)  =  120(1  -  ^  —}  ) 


(3.37) 


Ei  j  ydA  =  0 

yn  - 1 

yo  +  hi  yo+hi  +  hi 

j  EiydA  +      J  E2ydA  +        J  EiydA  =  0 


i  =  l 


yo  + J3i  +  h2  +  h3 


(3.38) 


yo 

where  hi 
h2 
h3 
El 

E2 


yo+hi 


yo  +  jji  +  h2 


the  thickness  of  first  layer (Steel) 
the  thickness  of  second  layer (Aluminum) 
the  thickness  of  third  layer (Steel) 
the  Modulus  of  material  1 (Steel) 
the  Modulus  of  material  2 (Aluminum) 


The  position  of  N.A.  (Yo)    can  be   calculated  by  using 
equation      (3.38).  Since      the     beam      is  symmetrically 

structured,  the  position  of  N.A.  must  be  at  the  center  of 
the  thickness  (Yo=  (hi+h2+h3) /2 )  .  The  flexural  rigidity  D  is 
calculated  as  follow: 
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iii^  hi  7 

D  =  t{Ei( —  +  hi(—  +  yo)  ] 
12  2 

h  ^  h 

+bLE2{—  +  h2{hi  +  —  +  yo)']  (3.39) 
12  2 

+Jb[Ei( —  +  h3(hi  +  h2  +  —  +  yo)  ] 
12  2 

where  b:  width  of  the  beam 

The  equivalent  moment  Mt(x)   is  obtained  as  follow: 

yo  +  hi  yo  +  hi  +  h2  yo  +  hi  +  h2  +  hi 


Mt  =    j  EiaiTydy  +     j  ETaiTydy  +        j  EiUiTydy  (3.40) 

yo  yo  +  iii  yo  +  hi  +  h2 

The  work  done  by  the  equivalent  moment  Mt(x)  is: 

L 

w  =  j  Mtv"dx 

0 

The  potential  energy  11  of  the  beam  and  the  constraints 
are  the  same  as  given  by  equations  (3.18)  and  (3.19).  The 
Lagrangian  function  L  is  also  the  same  as  expressed  in 
equation  (3.20).  Following  the  same  procedure  and  using 
equations  (3.21)  and  (3.22),  the  deflection  function  of  the 
beam  is  determined. 


>- 

>- 

>- 

> 

> 


Figure  3.2  6  The  F.E.A.  model  for  example  3.7 
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The  F.E.A.  model  used  for  comparison  is  shown  in  figure 
3.26.  This  model  contains  77  nodes  and  60  elements.  The 
results  from  the  proposed  method  for  different  number  of 
terms  N  are  shown  in  Figure  3.27.  The  assumed  deflection 
function  converges  at  N=7 .  The  comparison  between  the 
proposed  method  and  the  F.E.A.  is  shown  in  Figure  3.28.  The 
difference  in  the  maximum  deflection (x=27 )  is  approximately 
6.38%.  The  stress  distributions  through  the  thickness 
obtained  by  the  F.E.A.  and  proposed  method  are  shown  in 
figure  3.29.  The  maximum  shear  stress  at  the  interface 
surface  of  the  two  materials  as  obtained  from  the  proposed 
method  is  67.2  Mpa.  The  CPU  time  for  calculating  deflection 
by  the  F.E.A.  is  12  sec.  The  proposed  method  requires  less 
than  2  sec . 


Figure  3.27  the  computed  deflections  for  N=3 , 5 , 7 , 9 
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Figure  3.28  Deflections  by  the  proposed  method(N=7) 

and  the  F.E.A. 


The  results  illustrate  that  the  proposed  method  is  much 
faster  than  the  F.E.A.  and  that  the  two  methods  give  very 
close  deflection  values. 


7  (Tx 


-5.  10 


Figure  3.29  Stress  through  the  thickness  of  the  beam 
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Example  3 . 8 

Consider  a  clamped- supported  composite  beam  as  shown  in 
figure  3.30  subjected  to  the  temperature  distribution  T(x,y) 
as  shown  in  equation  (3.37).  The  composite  beam  in  this 
case  is  the  same  as  the  beam  of  example  3-7. 


rhi=i 

r-h2=l 

  L=30   * 

h3=l-' 

Figure  3.3  0  The  configuration  of  the  beam  in  example  3.8 

The  same  procedure  as  in  examples  3.7  is  followed  to 
obtain  Mt.  By  choosing  the  same  deflection  function,  the 
formulation  of  the  total  potential  energy  is  the  same  as 
that  in  example  3.2.  Since  the  boundary  conditions  are  not 
all  satisfied,  the  unsatisfied  boundary  condition  (v'=0  at 
x=0)  can  be  written  as: 

N 

G  =  ^  Ann  =  0 

n  =  l 

The  same  procedure  as  in  example  3.2  is  applied  to 
solve  for  the  An  and  X.  The  F.E.A.  model  used  for  comparison 
is  shown  in  Figure  3.31.  Since  the  boundary  conditions  and 
the  thermal  loading  are  not  symmetrical  in  the  x  direction. 


75 


this  model  contains  the  full  configuration  of  the  beam  with 
147  nodes  and  120  elements. 


Figure  3.31  The  F.E.A.  model  for  example  3.8 


The  results  from  the  proposed  method  for  different 
number  of  terms  N  are  shown  in  Figure  3.32.  The  assumed 
deflection  function  converges  at  N=10.  The  comparison 
between  the  proposed  method  and  the  F.E.A.  is  shown  in 
Figure  3.33.  The  difference  in  the  maximum  deflection (x=20 ) 
is  approximately  1.05%.  Since  the  temperature  loading  is 
the  same  as  in  example  3-7,  the  stress  distribution  through 
the  thickness  of  the  beam  and  shear  stress  is  the  same  as 
that  in  example  3.7.  The  CPU  time  for  calculating 
deflection  by  the  F.E.A.  is  11  sec.  The  proposed  method 
requires  less  than  1  sec. 
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5  10  15         20         25  30 


Figure  3.32  the  computed  deflections  for  N=5,10,20 


Again,  the  results  illustrate  that  the  proposed  method 
s  much  faster  than  the  F.E.A.  and  that  the  two  methods  give 
Imost  identical  deflection  values. 
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Example  3 . 9 

In  this  example,  the  composite  beam  has  the  same 
boundary  condition  as  in  example  3 . 6  and  the  same 
configuration  as  in  example  3.7.  The  beam  in  this  case  is 
shown  in  Figure  3.34. 


r-hl=l 

  L=3  0  

j-h2=l 

t  i 

h3=l-^ 

Figure  3.34  The  configuration  of  the  beam  in  example  3-9 


Following  the  same  procedure  as  in  examples  3-6  and  3- 
7,  the  deflection  is  obtained.  The  F.E.A.  model  used  is  the 
same  as  in  example  3 . 8  with  a  change  in  the  boundary 
conditions.  The  results  from  the  proposed  method  for 
different  niomber  of  terms  N  are  shown  in  figure  3.35.  The 
assumed  deflection  function  converges  at  N=8 .  The 
comparison  between  the  proposed  method  and  the  F.E.A.  is 
shown  in  Figure  3.36.  The  difference  in  the  maximum 
deflection (x=3  0)  is  approximately  1.49%.  The  CPU  time  for 
calculating  deflection  by  the  F.E.A.  is  12  sec.  The 
proposed  method  requires  less  than  2  sec. 
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Figure  3.35  the  computed  deflections  for  N=4,6,8 


Again,  the  results  illustrate  that  the  proposed  method 
is  much  faster  than  the  F.E.A.  and  that  the  two  methods  give 
almost  identical  deflection  values 


Figure  3.3  6  Deflections  by  the  proposed  method (N=10) 

and  the  F.E.A. 
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The  Two-Dimensional (Plate)  Problem 

In  this  investigation,  a  plate  is  defined  as  an  elastic 
body  bounded  by  two  flat  surfaces.  The  distance  between  the 
surfaces  is  the  thickness  of  the  plate.  The  neutral  surface 
(N.S.)  is  a  plane  parallel  to  the  top  and  bottom  surfaces. 
The  position  of  the  N.S.  can  be  obtained  in  the  same  fashion 
as  in  the  beam  case.  The  development  of  the  theory  of 
bending  of  a  thin  plate  is  based  on  the  following 
ass\imptions  attributed  to  Kirchhof  f  [20]  : 

1.  The  linear  filament  of  the  plate  which  is  initially 
normal  to  the  N.S.  remains  straight  and  normal  to 
N.S.  after  bending. 

2.  There  is  no  deformation  in  the  N.S. 

3 .  Normal  stress  in  the  direction  transverse  to  the 
N.S.  is  small  compared  to  the  other  normal 
stresses,  and  the  normal  strain  in  that  direction 
is  small  enough  to  be  disregarded. 

The  position  of  N.S.  for  a  composite  plate  shown  in 
figure  3.37  is  obtained  by  the  same  procedure  as  that  used 
in  the  beam  case.  Therefore,  the  total  flexural  rigidity  D 
of  the  composite  plate  is  defined  by  equation  (3.41).  The 
strain  energy  U  for  a  plate  is  formulated  as  equation 
(3 .42)  . 
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2  Ekik 
D  =  I 


(3.4i: 


where  Ik 
dk 

V 


:  hkdk^ 
1  , 

=    Vk  -  1  hk 

2 

:  Poisson's  ratio 

av    av.,         ,  ,av  av 


a' 


[7  =       I  {  (— -  +  —f  -  2(1  -  -U)  [—  —  -  (t^)  ]  }dA 


2     a'x  av 


a^x  a  y  dxdy 


(3.42) 


where  v  :  deflection  function 

A  :   the  surface  area  of  the  plate 


Due  to  the  effect  of  Poisson's  ratio,  the  equivalent 
moment  Mt  is  shown  in  equation  (3.43)  .  The  work  done  by  the 
equivalent  moment  Mt  is  given  in  equation  (3.44). 


h, 

 z 

 Zt 

h„ 

f 


N.S. 


Figure  3.37  The  cross  section  of  a  composite  plate 


Mt{x,  y)  =  ^  j  (1  +  "^^^  {EiaiT(x,  y,  z)z)dz 


(3.43) 


■i  =  l  Zi- 
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W  = 


(3.44) 


The    potential    energy    therefore    can    be    written  as 
follows : 


The  boundary  conditions  at  the  edges  of   the  plate  can 
be  categorized  as  follows: 

1.  clamped  edge  along  which  v  and  3v/3n  vanish 

2.  simply- supported  edge  along  which  v  and  Mn  vanish 

3 .  free  edge  along  which  Mn  and  Vn  vanish 

In  the  above  statements,  the  derivative  with  respect  to 
n  is  the  derivative  with  respect  to  the  normal  to  the  edge, 
and  Mn  and  Vn  denote,  respectively,  the  moment  and  the  shear 
force  in  the  normal  direction  to  the  edges.  The  procedure  of 
the  proposed  method  is  the  same  as  in  the  beam  case. 

Illustrative  examples 

Seven  examples  are  used  with  different  boundary 
conditions  in  order  to  demonstrate  the  capability  of  the 
proposed  method.  Examples  3.10  to  3.14  are  unit  homogeneous 
square  plate  with  dimension  Lx=Ly=l.  In  example  3.14,  a 
composite  plate  is  considered  in  order  to  demonstrate  both 
the  shear  stress  and  deflection.  In  examples  3.15  and  3.16, 
the   homogeneous    plates   with   holes    are    considered.    In  all 


n  =  U  -  W 


(3.45) 
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examples,  the  length  and  deflection  units  are  in  meters  and 
the  temperature  units  are  in  °C. 

In  order  to  illustrate  the  accuracy  and  the 
computational  efficiency  of  the  proposed  method,  it  is 
compared  with  the  finite  element  method  for  all  the 
considered  examples.  The  thermal  loading  for  these  seven 
cases  is  assumed  to  be  the  same.  The  temperature  difference 
between  top  and  bottom  surfaces  is  120°  and  the  temperature 
distribution  through  the  thickness  of  the  beam  is  assumed  to 
be  linear. 
Example  3.10 


s.s. 


s.s. 


Figure  3.3  8  the  configuration  of  the  plate  with  boundary 
conditions  in  example  3.10 


The  problem  of  a  homogeneous  plate  with  all  edges 
simply  supported  has  been  cited  by  many  authors  (for  example 
Timoshenko [20]  ,    Saada[27],    and  Love[28])    to  demonstrate  the 
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usage  of  Rayleigh-Ritz  method  in  plate  problems.  Taking  a 
unit  square  homogeneous  plate (thickness  h=0.01)  as  an 
example,  the  assumed  deflection  function (v)  is  chosen  to  be 
a  double  sine  series  for  both  x  and  y  direction. 

v{x,  y)  =  ^  ^  Am  sin  (mtx)  sin  (nny) )  (3.46) 

m=l n=l 

To  simplify  the  problem,  Nm  and  Nn  will  be  taken  the 
same  throughout  this  study.  With  the  deflection  surface 
approximated  by  the  double  sine  series  shown  in  equation 
(3.46),  the  strain  energy  U  which  is  obtained  by  equation 
(3.42)   is  expressed  as  : 

U  =  —  £  £A™„'(in'  +  n'f  (3.47) 

S     in  =  l  n  =  l 

The  temperature  distribution  T(z)  is  as  follow: 

nz)  =  120(1  -      "*"       )  (3.48) 
h 

The  equivalent  moment  Mt  and  The  work  done  by  the 
equivalent  moment  Mt  can  be  calculated  from  equations  (3.43) 
and  (3 .44) . 

120Da(l  -  1)') 

Mt{x,y)  =   (3.49) 

h 

V  V  ,  r..        +       (1  -  cos(in7t)  )  (1  -  cos(n7r)  )  ,^ 

^  =   \  \  Amn[Mt   (3.50) 

^inTi  rnn 

Since  the  second  derivative  of  a  sine  function  vanishes 
at  the  edges  of  the  plate,  i.e.  3V/9x^=0  at  x=0  and  x=l.  and 
3V/8y^=0    at    y=0    and    y=l.,     all    boundary    conditions  are 
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satisfied  by  the  assumed  deflection  function  and  no 
additional  constraint  equations  are  necessary.  The 
Lagrangian  function  L  is  equal  to  the  potential  energy 
function  in  this  case. 

Using  the  same  approach  as  in  the  beam  examples,  the 
partial  derivatives  of  L  with  respect  to  each  A^n  are  set 
equal  to  zero.  The  unknown  Amn  coefficients  are  obtained 
as : 

480a(l  +  D)  (1  -  cos  inm)  )  (1  -  cos  nK)  )  s 

Ama    -   -.   (i.bl) 

n  hmn 

In  view  of  the  above  equations,  it  can  be  seen  that  the 
numerical  value  of  ?^  diminishes  very  quickly  as  m  and  n 
increase . 


s.s. 


s.s. 


s.s. 


s.s. 


Figure  3.39  The  F.E.A.  model  for  example  3.10 


The  F.E.A.  model  used  for  comparison  is  shown  in  figure 
3.39.   This  model  contains  882  nodes  and  400  elements.  The 
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results  from  the  proposed  method  for  different  number  of 
terms  N  are  shown  in  figure  3.40.  Figure  3.41  shows  the 
comparison  between  the  F.E.A.  model  and  the  proposed  method 
for  the  deflection  curves  taken  along  the  section  y=0.5.  The 
assumed  deflection  function  converges  at  N=7 .  The 
difference  in  the  maximum  deflection  (x=0.5)  is 
approximately  0.56%.  The  deflection  of  the  plate  is  shown  in 
figure  3.42.  The  CPU  time  for  calculating  deflection  by  the 
F.E.A.  is  167  sec.  Th^e  proposed  method  requires  less  than  1 
sec . 


0.014 
G.G12 
0.01 
0.008 
0.006 
0.004 
0.002 

0.2  0.4  0.6  0.8  1 

Figure  3.40  the  computed  deflection  function  along  the 

section  y=0 . 5 
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V 


Figure  3.41  the  comparison  between  the  F.E.A.  model 
and  the  proposed  method (N=7)   for  the  deflection  curves 
along  the  section  y=0.5. 

Example  3 . 11 

As  illustrated  in  the  previous  section,  the  proposed 
method  makes  it  very  convenient  to  convert  one  plate  problem 
to  another  with  different  boundary  conditions.  It  is  also 
easy  to  change  the  boundary  conditions  of  a  plate  problem 
without  any  major  change  in  the  formulation  when  using  the 
proposed  method.  Consider  a  unit  square  homogeneous  plate 
with  three  edges  simply  supported  and  the  fourth  edge 
clamped  as  shown  in  figure  3.43. 
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y 


Figure  3.43  the  configuration  of  the  plate  with  boundary 

conditions  in  example  3.11 

With  the  same  assumed  deflection  function  as  equation 
(3.46),  the  formulations  of  the  strain  energy  and  the  work 
done  by  the  equivalent  moment  will  be  the  same  as  for 
example  3.10.  Since  the  boundary  conditions  are  not  all 
satisfied  by  the  assumed  deflection  function,  the  Lagrangian 
function  L  is  not  only  potential  energy.  The  unsatisfied 
boundary  condition  is  that  the  slop  at  the  fixed  boundary- 
must  be  zero. 

—  =  0         for  x=0  (3.52) 
ox 

Since   it   is   extremely  difficult  to   find  a  continuous 

function  to  force  the  slop  to  be  zero  along  the  entire  edge 
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x=0,  the  boundary  conditions  are  only  satisfied  using 
Lagrange  Multipliers  at  a  finite  number  of  points  chosen 
along  the  edge  x=0.  The  boundary  conditions  on  the  segments 
between  these  points  are  disregarded.  In  this  example,  the 
number  of  the  point  is  chosen  as  9  to  provide  enough 
accuracy  for  the  deflection  values. 

The  constraint  equations  are  obtained  as  follow: 

Gi(Am)  =  ^  ^  Amn(in7C)  sin(n7tyi)  =0   i=l  to  9  (3.53) 

m  =  l  m  =  l 

where  yi=i/10. 

The  Lagrangian  function  L  is: 

9 

L(Ann,  Xi)    =    U   -  W   +   ^    XiGilAmn)  )  '  (3  .  54) 

i  =  l 

Equating  the  partial  derivatives  of  L  with  respect  to 
each  Aon  to  zero  and  combining  the  result  with  equation 
(3.53),  there  are  (N^+9)  linear  equations  for  {N^+9) 
unknown  coefficients.  The  deflection  function  is  therefore 
determined. 

The  F.E.A.  model  used  for  comparison  is  the  same  as  in 
example  3.10  with  a  change  in  the  boundary  conditions.  The 
results  from  the  proposed  method  for  different  number  of 
terms  N  are  shown  in  figure  3.44.  Figure  3.45  shows  the 
comparison  between  the  F.E.A.  model  and  the  proposed  method 
for  the  deflection  curves  taken  along  the  section  y=0 . 5 . 
The    assumed   deflection    function    converges    at    N=17 .  The 
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difference  in  the  maximum  deflection  (x=0.65)  is 
approximately  0.90%.  The  deflection  of  the  entire  plate  is 
shovm  in  figure  3.46.  The  CPU  time  for  calculating 
deflection  by  the  F.E.A.  is  156  sec.  The  proposed  method 
requires  less  than  3  sec . 


Figure  3.44  the  computed  deflection  function  along 

the  section  y=0.5 
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Figure  3.45  the  comparison  between  the  F.E.A.  model  and 
proposed  method (N=17)   for  the  deflection  curves  along  the 

section  y=0 . 5 . 


Example  3 . 12 

In  this  example,  the  boundary  conditions  of  the  unit 
square  homogeneous  plate  are  2  simply  supported  edges  along 
x=0  and  y=0  and  two  clamped  edges  along  x=l  and  y=l.  The 
formulations  of  the  strain  energy  U  and  the  work  done  by  the 
equivalent  moment  will  be  the  same  as  in  previous  example. 
The  constraint  equations  are  obtained  as  follow: 
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N  N 


Gi{Aim)  =  ^  ^  Am{irm)  sin(nnyi)  =0,   1=1  to  9  (3.55a) 


N  N 


Gj(Amn)  =  ^  ^  Aminn)  sin  {nnxj  -  9)  =0,   j=10  to  18  (3.55b) 


jn~l  m  =  l 


where  Xi=(i/10),  yi=(i/10),  i=l  to  9 
The  Lagrangian  function  L  is : 

18 

UAnn.  'ki)    =   U   -  W   +   ^  XiGiiAm)  ) 

1=1 


(3.56) 


s.s. 


irlamped 


S.S, 


Clamped 


Figure  3.47  the  configuration  of  the  plate  with  boundary 
conditions  in  example  3-12 


Equating  the  partial  derivatives  of  L  with  respect  to 
each  Anm  to  zero  and  combining  with  equations  (3.55a)  and 
(3.55b),       there    are    (N^+18)    linear    equations    for  (N^+18) 


unknown     coefficients.  The     deflection     is  therefore 

determined.  .  ''•* 

The  F.E.A.  model  used  for  comparison  is  the  same  as  in 
example  3.10  with  a  change  in  the  boundary  conditions.  The 
results  from  the  proposed  method  for  different  number  of 
terms  N  are  shown  in  figure  3.48.  Figure  3.49  shows  the 
comparison  between  the  F.E.A.  model  and  the  proposed  method 
for  the  deflection  curves  along  the  section  y=0 . 5 .  The 
assumed  deflection  function  converges  at  N=17 .  The 
difference  in  the  maximum  deflection  (x=0.7)  is 
approximately  0.9%.  The  deflection  of  the  entire  plate  is 
shown  in  figure  3.50.  The  CPU  time  for  calculating 
deflection  by  the  F.E.A.  is  156  sec.  The  proposed  method 
requires  less  than  5  sec.  .  •  ■, 


V 


Figure  3.48  the  computed  deflection  function  taken  along 

the  section  y=0.5 
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Figure  3.49  the  comparison  between  the  F.E.A.  model 
and  the  proposed  method (N=17)   for  the  deflection  curves 

along  the  section  y=0.5. 

Example  3.13 

In  this  example,  the  boundary  conditions  of  the  unit 
square  homogeneous  plate  are  2  simply  supported  edges  along 
x=0  and  x=l  and  two  clamped  edges  along  y=0  and  y=l.  The 
formulations  of  the  strain  energy  U  and  the  work  done  by  the 
equivalent  moment  will  be  the  same  as  in  the  previous 
example.  The  constraint  equations  are  obtained  as  follows: 

Gi{Am)  -  5^  5]  Am(;nJU)  sin(n7Uyi)  =0,   1=1  to  9  (3.55a) 
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Gi  +  i(Amn)  =  ^  ^  Ann(mJi:)  sin  (mTC)  sin  (nTtyj)  =0  (3.55b) 


n=l m=l 


s.s. 


irlamped 


Clampled 


S.S, 


Figure  3.51  the  configuration  of  the  plate  with  the  boundary 

conditions  in  example  3.13 


The  Lagrangian  function  L  is: 

18 

UAm.  Xi)    =   U  -  W  +  XiGi(Ann)  )  (3.56) 

i  =  l 

Equating  the  partial  derivatives  of  L  with  respect  to 
each  Aran  to  zero  and  combining  with  equations  (3.55a)  and 
(3.55b),  there  are  (N^+18)  linear  equations  for  (N^+18) 
unknown  coefficients.  The  deflection  function  is  therefore 
determined. 
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The  F.E.A.  model  used  for  comparison  is  the  same  as  in 
example  3.10  with  a  change  in  the  boundary  conditions.  The 
results  from  the  proposed  method  for  different  number  of 
terms  N  are  shown  in  figure  3.52.  Figure  3.53  shows  the 
comparison  between  the  F.E.A.  model  and  the  proposed  method 
for  the  deflection  curves  along  the  section  y=0.5.  The 
assumed  deflection  function  converges  at  N=17.  The 
difference  in  the  maximum  deflection  (x=0.5)  is 
approximately  7%.  The  deflection  of  the  entire  plate  is 
shown  in  figure  3.50.  The  CPU  time  for  calculating 
deflection  by  the  F.E.A.  is  148  sec.  The  proposed  method 
requires  less  than  5  sec. 


Figure  3.52  the  computed  deflection  function  taken  along 

the  section  y=0.5 
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Figure  3.53  the  comparison  between  the  F.E.A.  model 
and  the  proposed  method (N=17)   for  the  deflection  curves 
taken  along  the  section  y=0.5. 

Example  3.14 

In  this  example,  a  unit  square  composite  plate  with  all 
edges  simply  supported  is  considered.  The  temperature  of 
top  surface  of  the  plate  is  120  degrees.  The  temperature  of 
the  bottom  surface  of  the  plate  is  zero.  The  temperature 
distribution  through  the  thickness  is  assumed  linear. 

The  neutral  surface  of  the  plate  is  obtained  by  using 
the  same  procedure  in  the  beam  problems.  The  total  flexural 
rigidity  D  of  the  composite  plate  is  therefore  determined  by 
equation   (3.41) . 
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In  this  example,  since  the  materials  are  symmetrical 
about  N.S.  The  N.S.  is  the  central  surface  throughout  the 
plate . 


s.s. 


s.s. 


Figure  3.55  The  plate  for  example  3.14 
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h3 

Figure  3.56  The  cross  section  of  the  plate  for  example  3.14 


The   configuration  and  cross   section  of   the  plate  are 
shovm  in  figure  3.56  and  the  thermal  loading  for  the  plate 
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is  as  in  equation  (3.57).  The  equivalent  moment  Mt  can  be 
determined  by  substituting  in  equation  (3.43). 

T(z)  =  120(1  -  ^^^^))  (3.57) 
h 

V      EiaiTix,  y,  z)z 

Mt(x,  y)  =  y   dz  (3.58) 

fr{/        1  -  \)i 

1-1  Zl-l 

By  assuming  the  deflection  function  to  be  the  same  as 
in  the  previous  examples,  given  in  equation  (3.46),  the 
strain  energy  U,  the  work  done  by  the  equivalent  moment  Mt 
and  potential  energy  are  therefore  formulated  as : 

U  =  —  £  £  A»>'(/n'  +  n'f  (3.59) 

8      m=l n=l 
Lx  Ly         -v  2  -\  2 

W  =   \  \  Mt{-^  +  ^)dxdy  (3  .  60) 

n  =  u  -  w 

since  the  second  derivative  of  a  sine  function  vanishes 
at  the  edges  of  the  plate,  i.e.  9V/3x^=0  at  x=0  and  x=l.  and 
9V/5y^=0  at  y=0  and  y=l.,  all  boundary  conditions  are 
satisfied  by  the  assumed  deflection  function  and  no 
additional  constraints  are  necessary.  The  Lagrangian 
function  L  is  equal  to  the  potential  energy  function. 

Using  the  same  approach  as  in  example  3.10,  the  partial 
derivatives  of  L  with  respect  to  each  A^n  set  equal  to  zero. 
The  unknown  A^n  coefficients  are  obtained  as: 


The  F.E.A.  model  used  for  comparison  contains  484  nodes 
and  300  elements.  The  results  from  the  proposed  method  for 
different  number  of  terms  N  are  shown  in  figure  3.57. 
Figure  3.58  shows  the  comparison  between  the  F.E.A.  model 
and  the  proposed  method  for  the  deflection  curves  along  the 
section  y=0.5.  The  assumed  deflection  function  converges  at 
N=7 .  The  difference  in  the  maximum  deflection  (x=0.5)  is 
approximately  8.62%.  The  shear  stress  distribution  through 
the  thickness  of  the  plate  obtained  by  The  F.E.A.  and  the 
proposed  method  are  shown  in  figure  3.58.  ! 


0.2         0.4         0.6        0.8  1 

Figure  3.57  the  assiomed  deflection  function  along  the 

section  y=0 . 5 
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V 
0.002 


Figure  3.58  the  comparison  between  the  F.E.A.  model 
and  the  proposed  method (N=7)   for  the  deflection  curves 
along  the  section  y=0.5. 


ax    10'  Pa 


Figure  3.59  shear  stress  distributions  through  the  thickness 
of  the  plate  obtained  by  The  F.E.A.  and  proposed  method 
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The  CPU  time  for  calculating  deflection  by  the  F.E.A. 
is  110  sec.  The  proposed  method  requires  less  than  1  sec. 
The  deflection  of  the  plate  is  shown  in  figure  3.60. 

Example  3.15 

A  plate  with  arbitrarily  shaped  holes  is  a  common 
structural  element,  but  is  extremely  difficult  to  handle  by 
the  classical  theory  of  plates.  The  finite  element  method, 
so  far,  is  the  approach  that  analysts  generally  use. 
However,  the  results  from  the  F.E.A.  is  not  a  continuous 
function.  It  can  only  produce  the  stress  and  deflection  at 
some  discrete  points (nodes) .  In  order  to  obtain  the 
deflection  values  by  F.E.A.  which  match  the  same  density  of 
the  results  from  the  proposed  method,  F.E.A.  would  require 
considerably  more  CPU  time  than  the  proposed  method. 

In  this  example,  the  unit  square  homogeneous  plate  with 
clamped  square  hole  is  considered.  The  edges  of  the  hole 
are  clamped  and  the  boundary  of  the  plate  is  simply 
supported. 

In  the  proposed  method,  a  pseudo  plate  of  the  same  size 
as  the  hole  can  be  fitted  into  the  hole,  with  its  boundary 
clamped.  The  pseudo  plate  can  then  be  integrated  with  the 
original  plate,  inside  which  the  hole  is  located,  and 
analyzed  as  a  whole  with  the  same  procedure  shown  in  example 
3.10.  The  pseudo  plate  inside  the  hole  can  be  then 
disregarded. 
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Figure  3.61  The  configuration  of  the  plate  for  example  3.15 


By  assuming  the  same  deflection  function  as  equation 
(3.46),  the  formulation  of  the  Lagrangian  function  L  remains 
the  same.  The  constraints  are  zero  deflection  and  zero 
normal  slope  along  the  edges  of  the  hole.  In  this  example, 
14  points  on  the  edges  of  the  hole  are  chosen  for  the 
constraints. 

The  F.E.A.  model  used  for  comparison  is  shown  in  figure 
3.62.  This  model  contains  1040  nodes  and  960  elements. 
Figure  3.63  shows  the  convergence (N=21)  of  the  proposed 
method  for  the  deflection  curves  along  the  section  y=0.375. 
Figure  3.64  shows  the  comparison  between  the  F.E.A.  model 
and  the  proposed  method  for  the  deflection  curves    along  the 
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section  y=0.375.  The  difference  in  the  maximiam  deflection 
(x=0.88)   is  approximately  5.50%. 


Figure  3.62  The  F.E.A.  model  for  example  3.15 


V 


Figure  3.63  the  computed  deflection  function  along  the 

section  y=0 . 375 


0.2  0.4  0.6  0.8  1 

Figure  3.64  the  comparison  between  the  F.E.A.  model 
and  proposed  method (N=21)   for  deflection  curves 
taken  along  the  section  y=0.375. 

The  deflection  of  the  entire  plate  is  shown  in  figure 
3.65.  The  CPU  time  for  calculating  the  deflection  by  the 
F.E.A.  is  223  sec.  The  proposed  method  requires  less  than  4 
sec. 

Example  3.16  ■  ■ 

An  illustration  of  the  analysis  of  a  plate  with  a 
square  clamped  hole  is  giving  in  the  previous  example.  It 
is  a  relatively  simple  problem  because  the  clamped  edges  can 
provide  the  necessary  forces  to  satisfy  the  boundary 
conditions  and  no  additional  forces  are  necessary.  For  the 
case  of  a  hole  with  free  edges,  a  different  technique  is 
necessary. 


Ill 


Clamped 


ClaiiiDed 


0.375 


Figure  3.66  the  configuration  of  the  plate  for  example  3.16 


Consider  a  unit  square  homogeneous  clamped  plate  with 
the  same  configuration  and  subjected  to  the  same  thermal 
loading  as  in  example  3.15.  To  begin  with,  a  clamped  plate 
without  a  hole  is  considered.  The  temperature  loading  is 
transformed  to  an  equivalent  moment  Mt.  The  deflection  under 
the  equivalent  moment  Mt  is  equated  to  zero  by  assuming  that 
the  equivalent  moment  Mt  is  less  than  Euler  loading (critical 
bucking  loading) .  For  the  case  where  the  equivalent  moment 
Mt  is  greater  than  Euler  loading (critical  bucking  loading), 
the  plate  can  buckle.  Due  to  the  nonlinear  behavior  of 
buckling,  the  proposed  method  and  the  linear  F.E.A.  are 
unable  to  determine  the  buckling  deflection  of  the  plate. 
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The  work(Wp)  done  the  equivalent  moment  Mt  when  a  plate 
is  without  a  hole  is  zero  by  assuming  there  is  no  buckling. 
The  work(Wh)  done  by  the  equivalent  moment  Mt  over  the 
region  where  the  hole  is  located  is  as  follows: 


0.625  0.625 

Wh 

0.375  0.375 


=    J     \  Mt{v")dxdy  (3.62) 


Wp  =  0  (3.63) 
The  total  work  done  is : 

W=Wp-Wh  (3.64) 
The  constraints  of  the  plate  are: 

Gi  =  v'(0,  yi)  =  0  for  i=l  to  9  (3.65a) 
where  yi=i/10 

Gi  +  9  =  v'(l,  yi)  =  0  for  i=l  to  9  (3.65b) 
where  yi=i/10 

Gi  *  IB  =  v'(xi,0)  =  0   for  i=l  to  9  (3.65c) 
where  Xi=i/10 

Gi  +  27  =  v'(xi,l)  =  0   for  1  =  1  to  9  (3.65d) 
where  Xi=i/10 

By  assuming  the  same  deflection  function  as  equation 
(3.46),  The  strain  energy  U  is  the  same  as  in  example  3.10. 
The  deflection  surface  of  a  clamped  plate  subjected  to  the 
thermal  loading  is  zero  and  the  moments  along  the  edges  of 
the  hole  can  be  calculated.  Since  the  moments  along  the 
edges  of  the  hole  that  is  free  must  be  zero,  external 
moments  are  superimposed  at  the  edge  of  the  hole  to  simulate 
the  effect  of  a  hole. 
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The  F.E.A.  model  used  for  comparison  is  the  same  as  in 
example  3.15  with  different  boundary  conditions.  Figure  3.67 
shows  the  convergence (N=16 )  of  the  proposed  method  for  the 
deflection  curve  along  the  section  y=0.375.  Figure  3.68 
shows  the  comparison  between  the  F.E.A.  model  and  the 
proposed  method  for  the  deflection  curves  along  the  section 
y=0.375.  The  difference  in  the  maximum  deflection  {x=0.5)  is 
approximately  1.26%.  The  deflection  of  the  plate  is  shown  in 
figure  3.69.  The  CPU  time  for  calculating  deflection  by  the 
F.E.A.  is  512  sec.  The  proposed  method  requires  less  than  6 
sec . 


Figure  3.67  the  computed  deflection  f unction (N=12, 16,20) 
taken  along    the  section  y=0.375 
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Proposed  method 
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Figure  3.68  the  comparison  between  the  F.E.A.  model 
and  proposed  method (N=16)   for  the  deflection  curves 
along  the  section  y=0.375. 


CHAPTER  4 
ILLUSTRATIVE  EXAMPLES 

Introduction 

In  chapters  2  and  3,  the  transient  temperature  field  of 
the  structure  is  obtained  by  the  simplified  finite 
difference  method.  The  thermal  deflection  of  a  structure 
subjected  to  an  assumed  thermal  loading (temperature  field) 
is  also  determined  by  using  the  proposed  method.  In  order 
to  solve  the  deflection  of  a  structure  under  heat  input,  it 
is  necessary  to  combine  the  two  schemes  discussed  in  the 
previous  chapters  in  order  to  evaluate  the  deflection  of  the 
structure.  In  other  words,  the  thermal  elastic  problem  is 
treated  in  two  steps: 

(1)  calculate  the  temperature  field  (transient  thermal 
problem) 

(2)  transform  the  temperature  field  into  mechanical 
loading  and  then  calculate  the  deflection  of  the 
structure 

The  transient  thermal  analysis  has  already  been  treated 
in  chapter  2.     After  the  thermal  loading  (temperature  field) 
is    evaluated,     the    deflection    of     the    structure    can  be 
determined  by  the  same  technique  discussed  in  chapter  3 . 


116 


117 


There  are  six  examples  in  this  chapter.  Examples  4.1 
and  4 . 2  are  beam  problems  and  examples  4.3  to  4.6  are  plate 
problems.  In  order  to  illustrate  the  accuracy  and 
computational  efficiency  of  the  results  from  the  proposed 
method,  it  is  necessary  to  compare  the  results  and  CPU  times 
with  the  F.E.A.  method. 

Example  4 . 1 

In  this  example,  a  homogeneous  multi-span  beam  shown  in 
figure  4.1  is  subjected  to  a  heat  input  (10  Watt)  at  two 
locations  (x=0.3  and  0.7)  on  the  top  surface  of  the  beam. 
The  temperature  value  at  x=0.2  to  0.4  and  x=0.6  to  0.8 
corresponding  to  the  full  penetration  depth  0 . 1  can  be 
directly  obtained  through  the  results  of  example  2.3.  The 
temperature  of  the  rest  of  the  beam  is  zero  during  this 
period. 


\ 

^ — -1  — ^ 
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h=0.1 

r  0.3 
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Figure  4.1  The  configuration  of  the  beam  for  example  4.1 
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After  obtaining  the  temperature  f ield (T (x, y) )  of  the 
beam  at  that  time,  the  equivalent  thermal  loading (Mt)  can 
be  obtained  as : 

n  h/2 

Mt{x)  =  X  ^  J        y)y(^y  ( 4 .  i) 

J=l  -h/2 

Where    b:   the  width  of  the  beam 

The  strain  energy  can  be  expressed  as  function  of 
deflection (v)  as  given  in  equation  (4.2).  The  work  done  by 
the  equivalent  loading  (Mt)  also  can  be  presented  as  a 
function  of  deflection. 

f  EI 

U  =  \   {v"fdx  (4.2) 

1 

w  =  j  Mt(v")dx  (4.3) 

0 

The    boundary    conditions     of     the    beam    are    v=0  at 

x=0,0.6L,    L   and  3^v/3^x=0    at   x=0    and   L.       By   assuming  the 

deflection  function  v  as  a  Fourier  series: 

,rmx^               .    ,  nnx^ 
V  =  Ao  +  2^  An  cos  (  )  +  2^  Bn  sm  (  ) ) 

0=1  ^  n=l  ^ 

The   strain   energy  U   is    the   same   as   equation  (3.26). 
The  work  done  by  the  equivalent  loading (Mt  (x)  )  is: 


w 

n=l  "0  Ji=l 


=        Ann^n^j  Mt{x)  COS  {mix)dx  -  ^  Bnri^n^j  Mt{x)  s in  (nKx)dx 

0 

(4.4) 


The  constraints  are  thus: 


Gi  =  Ao  +  2,  An  =  0  (4 .  5a) 


n  =  l 


Gz  =  ^n^Axi  =  0  {4.5b) 


n  =  l 


G3  =  Ao  +  ^AnCos(n7l)  =  0  (4.5c) 


n=l 


G4  =  2]  A^n^  cos  (nTT)  =  0  (4.5d) 


"  .nTCXs.       ^  nnxs. 


An  COS  (  )  +  +2^BnCOS(  )  =  0  (4,5e) 


n=l  ^  n=l  ^ 


Where  Xs=0.6L 


The  Lagrangian (merit )  function  L  can  be  formulated  as: 


5 


L  =  U  -  w  +  ^  XiGi 


i  =  l 


By  evaluating  to  zero  the  partial  derivatives  of  the 
merit  function  with  respect  to  parameters  Xi  (1=1,2,3,4,5) 
Ao,  An  and  Bn  (n=l , 2 , . . .N) ,  a  set  of  (2N+6)  linear  equations 
with  (2N+6)  unknowns  is  obtained.  Therefore  all  the  unknown 
coefficients  can  be  solved  and  the  deflection  function  of 
the  beam  is  determined. 
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The  F.E.A.  model  used  for  comparison  contains  the  full 
configuration  of  the  beam  with  11111  nodes  and  1000 
elements.  The  results  from  the  proposed  method  for  different 
number  of  terms  N  are  shown  in  Figure  4.2.  The  assumed 
deflection  function  converges  at  N=10 .  The  comparison 
between  the  proposed  method  and  the  F.E.A.  is  shown  in 
Figure  4.3.  The  difference  in  the  maximiam  deflection  (x=0  . 3  ) 
is  about  7.47%.  The  CPU  time  for  calculating  deflection  by 
the  F.E.A.  is  7  6  seconds.  The  proposed  method  requires  less 
than  2  sec . 


Figure  4.2  the  computed  deflection  for  N=8,10,12 
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Figure  4.3  comparison  of  the  deflection  results  from  the 
proposed  method (N=10)  and  the  F.E.A. 

Again,  the  results  illustrate  that  the  proposed  method 
is  much  faster  than  the  F.E.A.  and  that  the  two  methods  give 
very  similar  deflections  values. 


The  composite  beam  shown  in  figure  4.4  is  assumed  to  be 
subjected  to  heat  flux  at  the  top  surface  of  the  beam.  The 
temperature  distribution  through  the  thickness  of  the  beam 
can  be  directly  obtained  from  the  results  of  example  2.2. 
The  temperature  through  the  depth  of  the  beam  is  independent 
of  X. 


Example  4 . 2 
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Figure  4 . 4  The  configuration  of  the  beam  for  example  4 . 2 


The  beam  is  clamped  at  x=0  and  simply  supported  at 
x=1.4.  The  position  of  N.A.  (yo)  and  the  total  flexural 
rigidity  D  can  be  obtained  by  equations  (3.33)  and  (3.34). 
The  equivalent  moment  Mt(x)   is  obtained  as  follows: 


yo  +  hi  yo  +  hi  +  h2 

Mt  =     j  EiUiTydy  +    j  EM2Tydy  (4.6) 

yo  yo  +  hi 


The  work  done  the  equivalent  moment  Mt(x)  is: 


w  =  j  Mtv"dx  (4.7) 


The  constraints  in  this  case  are: 


Gi  =  Ao  +  ^An  =  0  (4.8a) 

n  =  l 


G2  =  ^nBn  =  0  (4.8b) 

n  =  l 
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(4.8c) 


n  =  l 


N 


(4.8d) 


n=l 


Gs  =  Ao  +  2^  A,  COS  (  )  +  +2^  Bn  cos  (  )  =  0        (4  .  Be) 


Where  Xs=0.7L 

By  assuming  the  same  deflection  function  as  in  example 
4.1,  The  Lagrangian (merit)   function  L  can  be  formulated  as: 


By  equating  to  zero  the  partial  derivatives  of  the 
merit  function  with  respect  to  parameters  X,i  (1=1,2,3,4,5) 
Ao,  An  and  Bn  (n=l ,  2 ,  .  .  .N)  ,  a  set  of  (2N+6)  linear  equations 
with  (2N+6)  unknowns  is  obtained.  Therefore  all  the  unknown 
coefficients  can  be  evaluated  and  the  deflection  function  of 
the  beam  is  determined. 

The  F.E.A.  model  used  for  comparison  contains  the  full 
configuration  of  the  beam  with  441  nodes  and  400  elements. 
The  results  from  the  proposed  method  for  different  number  of 
terms  N  are  shown  in  Figure  4.5.  The  assumed  deflection 
function  converges  at  N=6 .  The  comparison  between  the 
proposed  method  and  the  F.E.A.    is  shown  in  Figure  4.6a.  The 


5 


i  =  l 


difference  in  the  maximiim  deflection (x=2 )  is  3.83%.  The  CPU 
time  for  calculating  the  deflection  by  the  F.E.A.  is  37 
seconds.     The  proposed  method  only  requires  less  than  1  sec. 


Figure  4.5  the  calculated  deflection  of  the  beam  for  N=4,6,8 

The  stress  distribution  through  the  thickness  also  can 
be  determined  by  the  following  procedure: 

Oni)  =  -0LiEiT  For  layer  i(i=l  to  2)  '       '  '. 

Caved)  =  —  ^  jaiEiTdy  For  layer  i(i=l  to  2) 

Eiv  7 

Omii)  =         >      CLjEjTydy  For  layer  i{i=l  to  2) 
D  fr[  ^ 

The  resultant  stress  for  any  layer  of  the  beam  can  be 
calculated  as  (5  =(5i+Guv,+Cm  .  The  stress  results  of  the 
proposed  method  and  the  F.E.A.  are  shown  in  figure  4.6b. 
The  slip(shear)   stress  as  shown  in  figure  4.6b  which  is  the 
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resultant  stress  difference  between  the  two  material  is  43 
Mpa . 


V 


-2.00Ef08   

thickness 


(b) 

Figure  4.6   (a)   deflection  by  the  proposed  method (n=6)  and 
the  F.E.A.    (b)  Resultant  stresses  a  through  the  thickness  of 

the  beam 


i 
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Again,  the  results  illustrate  that  the  proposed  method 
is  much  faster  than  the  F.E.A.  and  that  the  two  methods  give 
almost  identical  numerical  values . 


Example  4 . 3 

The  problem  of  a  homogeneous  plate  subjected  an  assiamed 
thermal  loading  with  all  edges  simply  supported  has  been 
solved  in  example  3.10.  A  unit  homogeneous  square  plate  as 
shown  in  figure  4.7  (thickness (h) =0 . 1)  is  considered  and 
subjected  to  a  heat  input (20  Watt)  at  x=0.3  y=0.3  of  the 
upper  surface.  The  assiamed  deflection  function (v)  is  chosen 
to  be  a  double  sine  series  for  both  the  x  and  y  direction. 


s.s. 


s.s. 


Figure  4.7  the  configuration  of  the  plate  with  boundary 
conditions  in  example  4.3 
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vix,  y)  =  X  X  ^        {ntltx)  sin  {rmy) 

m=l  n=l 

The    strain    energy    U    which    is    obtained    by  equation 

(3.42)  is  expressed  as  : 

S      jn  =  l  n  =  l 

The  temperature  distribution  T(x,y,z)  is  developed  as 
in  example  2.5.  The  equivalent  moment  Mt  and  The  work  done 
by  the  equivalent  moment  Mt  can  be  calculated  from  equations 

(3 .43)  and  (3 .44)  . 

0.05 

MtOc,  y)  =    J  EaT{x,  y,  z)  (1  +  v)zdz  (4.9) 

-0.05 

W  =  Mti—j  +  —)dA  (4.10) 

jii=i  n=i    OA  y 

since  the  second  derivative  of  a  sine  function  vanishes 
at  the  edges  of  the  plate,  i.e.  3V/8x^=0  at  x=0  and  x=l.  and 
3^v/9y^=0  at  y=0  and  y=l.,  all  boundary  conditions  are 
satisfied  by  the  assumed  deflection  function  and  no 
additional  constraint  equations  are  necessary.  The 
Lagrangian  function  L  is  equal  to  the  potential  energy 
function  in  this  case. 

Using  the  same  approach  as  in  the  beam  examples,  the 
partial  derivatives  of  L  with  respect  to  each  A^n  are  set 
equal  to  zero.  The  unknown  A^n  coefficients  are  obtained 
as : 
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4f   Wt{x,  y)  sin  {mix)  sin  (miy)dA 
Amn  =   (4.11) 


K^Ddii^  +  n^) 


Due  to  the  limitation  of  the  number (32000 )  of  the  nodes 
in  F.E.A.  software,  it  is  not  possible  to  evaluate  the 
temperature  in  the  plate  with  equally  space  mesh.  In  order 
to  predict  a  converge  solution,  three  different  F.E.A. 
models  are  used  with  the  following  structure: 

1.  20  elements  in  the  x  direction,  20  elements  in  the  y 
direction  and  2  elements  in  the  z  direction. 

2.  40  elements  in  the  x  direction,   40  elements  in  the  y 
direction  and  4  elements  in  the  z  direction. 

3.  50  elements  in  the  x  direction,   50  elements  in  the  y 
direction  and  5  elements  in  the  z  direction. 

The  maximum  deflections  (at  x=0.3,  y=0.3)  and  CPU 
time  for  the  three  different  F.E.A.  models  and  the  results 
from  the  proposed  method  are  shown  in  table  4.1.  The 
estimated  converging  value  from  the  F.E.A.  models  using  a 
very  large  number  of  elements  for  the  maximum  deflection  is 
approximately  the  same  as  that  calculated  by  the  proposed 
method  (see  figure  4.8a), 


'F.E.A.  CPU  time 
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1                 Proposed  method: 

deflection:  9.7x10"^ 

CPU  time:  2  sec 
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(b) 


Figure  4.8   (a)   the  comparison  of  the  results  from  the 
proposed  method  and  the  F.E.A.    (b)   the  calculated  deflection 
function  along    the  section  y=0,3 
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Table  4.1 


the  CPU  time 

the 

node  / 

( sec ) 

deflection  at 

element  # 

x=0.3 

y=0.3 

F.E.A.  model  1 

386 

16.4x10"' 

1323/800 

F.E.A.  model  2 

2930 

14.12x10'' 

8405/6400 

F.E.A.  model  3 

11588 

13  .23x10'' 

15606/12500 

The  proposed 

2 

9.7x10"' 

method (N= 9) 

The  results  from  the  proposed  method  and  the 
extrapolated  F.E.A.  results  (figure  4.8a)  are  almost 
identical.  The  results  from  the  proposed  method  with 
different  niimber  of  terms  N  are  shown  in  figure  4.8b.  The 
assumed  deflection  function  converges  at  N=ll.  The  CPU  times 
for  calculating  deflection  by  the  F.E.A.  are  orders  of 
magnitude  larger  than  those  for  proposed  method. 


Example  4 . 4 

In  this  example,  the  homogeneous  unit  square  plate  with 
clamped  square  hole  is  considered.  The  edges  of  the  hole 
are  clamped  and  the  boundary  of  the  plate  is  simply 
supported.  The  plate  with  boundary  conditions  and  subjected 
to   the   same  heat   input   in  example   4.3    is   shown   in  figure 
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Following  the  same  procedure  as  in  example  3.15,  a 
pseudo  plate  of  the  same  size  as  the  hole  can  be  fitted  into 
the  hole,  with  its  boundary  clamped.  The  pseudo  plate  can 
then  be  integrated  with  the  original  plate,  inside  which  the 
hole  is  located,  and  analyzed  as  a  whole  with  the  same 
procedure  as  in  example  3.15.  The  pseudo  plate  inside  the 
hole  can  be  then  disregarded. 


s.s. 


0.2 

r 

s.s. 

Clanped 

Qin 

0.2 

J 

1 

0.4 

1 

S.S. 

0.4 

Figure  4.9  The  configuration  of  the  plate  for  example  4.4 

By  taking  the  same  deflection  function  as  in  equation 
(3.46),  the  formulation  of  the  Lagrangian  function  L  is  the 
same  as  in  example  4.3.  The  constraints  are  zero  deflection 
and  zero  normal  slope  along  the  edges  of  the  hole.     In  this 
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example,  14  points  are  chosen  on  the  edges  of  the  hole  for 
imposing  the  constraints. 

Figure  4.10  shows  the  convergence (N=21)  of  the  proposed 
method  for  the  deflection  curves  along  the  section  y=0.3. 


Figure  4.10  the  calculated  deflection  function  along  the 

section  y=0 . 3 


Due  to  the  accuracy  of  the  results  from  the  proposed 
method  as  demonstrated  in  example  3 . 15  and  checked  by  the 
F.E.A.  model  and  the  extensive  computations  necessary  for 
the  F.E.A.  model  to  simulate  the  exact  temperature  field  of 
the  plate,  the  F.E.A.  model  is  omitted  in  the  rest  of  the 
examples  in  this  chapter.  The  convergence  of  the  proposed 
method  will  be  considered  as   the   solution  without  further 
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checking.  The  proposed  method  requires  less  than  5  sec  for 
obtaining  the  converging  solution. 

Example  4 . 5 

In  this  example,  a  unit  square  composite  plate  with 
mixed  boundary  condition  is  considered.  Since  the  proposed 
method  allows  the  boundary  conditions  to  be  satisfied 
through  the  use  of  Lagrange  Multipliers  on  a  finite  number 
of  discrete  points  along  the  boundary,  it  is  not  difficult 
to  handle  a  mixed  boundary  condition  problem  without  any 
major  modification  of  the  assumed  deflection  or  energy 
functions . 

The  strategy  of  applying  the  proposed  method  to  a  mixed 
boundary  condition  problem  is  to  impose  suitable  constraints 
according  to  the  desired  boundary  conditions  at  selected 
discrete  points  on  the  boundaries  which  are  not  completely 
satisfied  by  the  assumed  deflection  function.  The  plate 
with  the  boundary  conditions  and  thermal  loading  is  shown  in 
figure  4.11.  The  temperature  field  is  the  same  as  obtained 
in  example  4.2. 

By  assuming  the  same  deflection  function  and  applying 
the  same  procedure  as  in  the  previous  examples,  the  strain 
energy  U,  the  work  done  by  the  equivalent  moment  Mt  and  the 
potential  energy  are  therefore  determined.  The  additional 
constraint  points  are  chosen  as  the  following  sets  of  points 
along  the  edges  of  the  plate: 
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set  1:  3v/3x=0  at  x=0  and  y=0.1, 0.2, 0.3 
set  2:  9v/3x=0  at  x=1.0  and  y=0 . 1 , 0 . 2 , 0  .  3 


set  3:  3v/dy=0  at  y=0  and  x=0 . 1 , 0 . 2 , 0  .  3 


set  4:  3v/3y=0  at  y=1.0  and  x=0 . 1 , 0 . 2 , 0  .  3 


Figure  4.11(a)   the  plate  with  boundary  condition  for  example 
4.5   (b)   the  thermal  loading  plate  for  example  4.5 
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Therefore,  The  Lagrangian  function  L  is  the  combination 
of  the  potential  energy  function  and  14  constraints.  Using 
the  same  approach  as  in  example  4.4,  The  unknown  A^n 
coefficients  are  obtained.  The  results  from  the  proposed 
method  for  different  niomber  of  terms  N  are  shown  in  figure 
4.12a.  The  assumed  deflection  function  converges  when  N=19 
to  21.  The  proposed  method  requires  less  than  5  sec  to 
obtain  a  converging  solution. 

The  resultant  stress  for  any  layer  of  the  beam 
direction  can  be  calculated  as: 

O  =Ol+Oave+<5m 

where 

Gt(i)  =  -(1  +  VilCLiEiT  for  i  =  1  to  2 

Caved)  =  — — - — ^  JaEiTd  +  Vi)dz  for  i  =  1  to  2,  n=2 

X  hjEj  — 

EiV  ^  7 

Gmii)  =         >      ajBjTd  +  Mi)zdz  for  i  =  1  to  2 ,  n=2 

J-l  ZJ-l 

Because  the  thermal  loading  is  the  same  as  x  and  y 
direction,  the  stresses  in  x  and  y  directions  (ax,ay)  are 
the  same  as  the  resultant  stress.  The  shear  stress (Op)  can 
determined  by  combined  the  ax  and  Qy. 
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(b) 

Figure  4.12   (a)   the  calculated  deflection  function  along 
the  section  y=0.5   (b)  Resultant  stress  through  the  thickness 

of  the  plate 
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The  stress  results  from  the  proposed  method  and  the 
F.E.A.  are  shown  in  figure  4.12b.  The  slip(shear)  stress  as 
shown  in  figure  4.12b  which  is  the  difference  of  the 
resultant  stress  between  the  two  material  is  84.5  Mpa. 

Example  4 . 6 

An  illustration  of  how  to  analyze  a  plate  with  a  square 
free  hole  is  given  in  example  3.16.  Consider  a  clamped 
composite  plate  with  a  central  square  hole  to  be  subjected 
to  the  same  thermal  loading  given  in  example  4.5  as  shown  in 
figure  4.13.  The  cross  section  of  the  plate  is  also  assumed 
to  be  the  same  as  in  example  4.5. 


Clamped 


ClamDed 


0.375 


Figure  4.13  the  configuration  of  the  composite  plate  for 

example  4 . 6 
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Following  the  same  procedure  as  in  examples  4 . 5  and 
3.16,  the  temperature  loading  is  transformed  to  an 
equivalent  moment  Mt(x,y)  and  the  work(Wp)  done  by  the 
equivalent  moment  Mt  and  the  work(Wh)  done  by  the  equivalent 
moment  Mt  over  the  hole  region  are  calculated. 

0.62  5  0.62  5 

Wh  =    j     J  Aft(7C,  y)v"dxdy  (4.12) 

0.375  0.375 

Wp  =  0  (4.13) 
The  total  work  done  is : 

W=Wp-Wh  (4.14) 
The  rest  of  the  procedure  of  obtaining  the  deflection 
function  is  the  same  as  in  example  3.16.  Figure  4.14  shows 
the  convergence (N=16 )  of  the  proposed  method  for  the 
deflection  curve  along  the  section  y=0.375.  The  CPU  time 
for  calculating  the  deflection  function  by  the  proposed 
method  in  order  to  obtain  the  converging  solution  is  less 
than  5  sec . 
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CHAPTER  5 
CONCLUSIONS  AND  RECOMMENDATIONS 

Conclusions 

It  can  be  concluded  from  this  study  that  the  proposed 
method  provides  an  excellent  tool  for  the  analysis  of  a  wide 
variety  of  thin  elastic  structures  subjected  to  thermal 
loading. 

The  proposed  method  can  be  used  to  treat  almost  all 
type  of  thin  elastic  structures  with  different  heat  inputs 
and  boundary  conditions.  The  main  advantage  of  the  proposed 
method  is  that  the  boundary  conditions  can  be  changed 
without  major  modification  of  the  formulation  as  illustrated 
in  the  different  cases  considered  in  this  studies. 

The  proposed  method  also  can  be  applied  to  many  plate 
problems  with  discrete  and  mixed  boundary  conditions  as 
illustrated  in  example  4.5  and  for  plates  with  different 
shapes  and  types  of  holes.  In  all  these  cases,  the  boundary 
conditions  can  be  changed  without  major  modification  of  the 
formulation. 

Whenever  applicable,  the  proposed  method  has  several 
advantages  as  summarized  in  the  following: 
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The  proposed  method  provides  continuous  (closed 
form)  expressions  for  the  deflection  and  stress 
solution  while  the  F.E.A.  only  provides  discrete 
results  at  the  nodes  in  the  considered  mesh. 
The  F.E.A.  can  not  directly  provide  the  slip (shear) 
stress  at  the  common  surface  of  two  different 
materials  for  a  composite  material  structure  which 
is  subjected  to  a  thermal  loading.  The  proposed 
method  can  readily  determine  this  information.  ■  ^  - 
The  proposed  method  can  evaluate  the  converging 
solution  in  significantly  shorter  CPU  time  for  both 
thermal  loading  and  elastic  structural  analysis. 
Very  large  F.E.A.  models  with  fine  meshes  or 
special  mesh  distributions  may  be  needed  for  some 
problems  in  order  to  obtain  the  same  converging 
results  as  the  proposed  method.  For  example,  in 
example  4.3,  the  12500  equal  size  elements  model 
requires  about  3  hours  CPU  running  time  to  evaluate 
the  results  and  this  model  still  does  not  converge. 
On  the  other  hand,  less  than  5  seconds  are  required 
to  obtain  the  converging  solution  by  the  proposed 
method.  •  ■'. 
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Recommendations 
The    proposed    method,     although    possessing  valuable 
advantages    in   analyzing    thermal    deflection   and    stress  in 
thin     elastic     structural     elements,      has     the  following 
disadvantages : 

1.  The  finite  difference  segment  of  the  proposed 
method  can  only  obtain  numerical  temperature 
distributions  (not  continuous  functions) ,  the  size 
of  the  control  volume  in  the  transient  thermal 
analysis  becomes  very  important  in  providing  the 
closest  thermal  equivalent  loading. 

2.  For  a  clamped  beam  or  a  clamped  plate,  there  are 
zero  deflection  results  due  to  the  nature  of  the 
thermal  equivalent  moment.  However,  if  the  load  is 
close  to  or  greater  than  the  Euler  load,  the 
structure  undergoes  buckling  and  the  proposed 
method  cannot  deal  with  this  non-linear  condition 
to  determine  the  deflection. 

The  following  studies  are  recommended  for  future 
extension  of  the  method  presented  here: 

1.  It  is  suggested  that  an  efficient  procedure  should 
be  developed  in  order  to  apply  the  proposed  method 
to  problems  with  combined  thermal  loading  and 
mechanical    loading.      This    should  not   present  any 
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difficulties  since  superposition  can  be  used  in 
this  case. 

It  is  suggested  that  the  proposed  method  is 
extended  to  deal  with  non- linear  problems  such  as 
beams  or  plates  with  large  deflections  by  iterative 
or  piece-wise  linear  analysis  or  averaging 
techniques . 

It  is  suggested  that  a  procedure  is  developed  to 
apply  the  proposed  method  for  structures  where 
buckling  can  occur. 

It  is  suggested  to  extend  the  proposed  method  for 
treating  structures  subjected  to  fluctuing  thermal 
loading  in  order  to  predict  their  fatigue  life. 


APPENDIX 


DETAILED  ILLUSTRATION  OF  THE  PROCESS  OF  DETERMINING  THE 
THERMAL  DEFLECTION  AND  STRESS  FOR  A  COMPOSITE  BEAM 
SUBJECTED  TO  A  GIVEN  HEAT  INPUT 


Consider  a  composite  beam  with  boundary  conditions  and 
subjected  to  heat  input  (Qin=10  Watt  for  120  seconds)  at 
x=0.3  on  the  top  surface  of  the  beam  as  shown  in  figure  A.l. 
The  transient  temperature  distribution  for  the  beam  at  the 
end  of  the  heating  time  is  first  obtained  before  calculating 
the  thermal  deflection. 


( 

r 

  1  m  

Steel 

Aluminiom 

 0.3m  

A  J 

0.1  m 

 0.7m  

Figure  A.l  The  configuration  of  the  beam 


Width  of  the  beam,  b=0.01  m 
Material  1   :  Steel 

Density  of  material  ,pi=7900  (Kg/m^) 

Thermal  conductivity    Ki  =  45  (W/m-°K) 

Specific  heat  ci=460  (J/kg-°K) 

Young's  Modulus  of  elasticity  Ei=2 . 1x10^  (Pa) 

Thermal  diffusivity  pi=  Ki/ (piCi)  =1 . 238x10'^  (W-mV  J) 
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Material  2  :  Aliuninum 

Density  of    material , p=2700  (Kg/m^) 

Thermal  conductivity,   K  =  200  (W/m-°K) 

Specific  heat,   c=900  (J/kg-°K) 

Young's  Modulus  of  elasticity  E=7 . 9x10^  (Pa) 

Thermal  diffusivity  p2=  K2/ (P2C2)  =8  .  230x10"^  (W-mV  J) 

The  Temperature  Distribution 

The  temperature  distribution  for  the  beam  is  calculated 
by  using  the  simplified  finite  difference  method.  The  first 
step   is   to   calculate   the   surface   area  and  voliime   of  each 

control  volume  (layer)  by  selecting  the  thickness (Ah)  of 
the  incremental  layer  in  the  first  material  (Steel)  as  0.01 
m  to  allow  for  10  penetration  steps.  In  order  to  simulate 
the  temperature  at  the  position  of  the  heat  input,  the 
radius  of  the  control  volume  (layer)  should  has  much  smaller 
volume  than  the  other  layers  and  is  taken  as  half  of  Ah 
and  the  finial  temperature  of  layer  1  can  be  assumed  to  be 
the  temperature  at  the  position  of  the  heat  input. 

The  radius  and  required  heat  penetration  times 
(calculated  from  equation  (A.l))  of  the  control  volumes  are 
shown  in  table  A.l, 


n  ^ 

=  0.0891— 


(A.l) 
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Table 

A 

.1 

Control 
volume  # 

Dp 

tp 

1 

0. 

005 

0.1578 

2 

0. 

015 

1.4202 

3 

0. 

025 

3 .9449 

4 

0. 

035 

7 .7320 

5 

0  . 

045 

12 .7815 

6 

0  . 

055 

19.0933 

7 

0. 

065 

26.6675 

8 

0. 

075 

35.5041 

9 

0. 

085 

45.6031 

10 

0. 

095 

56 .9644 

11 

0. 

105 

69.5880 

12 

0. 

115 

93 .4741 

13 

0. 

125 

98.6225 

14 

0. 

135 

115.0333 

15 

0. 

145 

132.7065 

Qin 


Figure  A. 2  The  control  volumes  of  the  beam 
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The  surface  area  and  volume  of  each  control  volume  are 
shown  in  figure  A. 2.     They  are  calculated  as  following: 

The  voliame  and  surface  area  of  material   1   in  layer  1 

is : 

Vi,i  =    (A.  2) 

8 

bAhK 

Ai.i  =    (A.  3) 

2 

From  layers  2  to  10,    the  shape  of  the  surface  area  is 

half -cylindrical,    the  thickness  of  each  layer  is  Ah (=0.01) 

and  their  volxome  and  surface  in  material  1  are: 

jb(4j  -  DAh^n   ^       .  ^        ,  ^  ,     , , 

1^1,  J  =   for  j=2  to  10  (A.  4) 

2 

Ai,,  =  "  ^^^^    for  j=2  to  10  (A.  5) 

2 

Since  the  first  10  layers  do  not  involve  material  2, 
the  volumes  and  surface  areas  in  material  2  of  the  first  10 
layers  are  zero . 

V2.  j  =  0.  for  j=l  to  10  (A.  6) 

A2,  j  =  0.  for  j=l  to  10  (A. 7) 

The  Layer  11  as  shown  in  figure  A. 3  has  both  materials 

1    and    2 .       The    volume    and    surface    area    of    layer    11  in 

material  1  is: 


Vi.ii  =  Jb|^area_  oab  +  area_  ode  +  Aobc  -  ^ 


(A. 8) 
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Ai,  j  =  2JbDpii0ii 


(A. 9) 


Dpii  =  0.105 
DPin  =  0  •  095 


Figure  A. 3  The  layer  11 


Since  the  point  p  (in  material  2)  and  b  (in  material  1) 
are  on  the  same  heat  exchange  surface,  they  have  the  same 
penetration  distance  Dpn. 

Dpii  =  ri  +  r2  1-^  (A.  10) 

where  ri :  line_og=hiCsc  (9) 
r2:  line_gp 

From  equation  (A. 10),  the  function  of  curve_bpc  can  be 
determined  by  equation  (A. 11)  in  cylindrical  coordinate 
system: 

r  =  ri  +  r2  =  P  -  0  csc(e))  (A.  11a) 

where 
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P  =  Jgqpn 
0  =  (1  -  J|)h> 

It-en 

V2,ii  =  jb(  jr'de  -  AoJbc))  (A. lib) 

Su 

It  -Bu 

>i  A2,ii  =  curve_  bpc  =  b  j  rd0  (A.  11c) 

en 

The  volume  and  surface  area  in  material  2  of  layer  11 
can  be  calculated  through  the  integration  of  equations 
(A. lib)  and  (A. 11c).  The  volumes  and  surface  areas  in  the 
rest  of  the  layers  are  calculated  in  a  similar  way.  The 
volume  and  surface  area  for  each  layer  are  shown  in  table 
A.2 . 


Table  A.2 


control 

volume  in 

volume  in 

area  in 

area  in 

volume # 

material  1 

material  2 

material  1 

material  2 

1 

0.0001963 

0 

0.00157 

0 

2 

0.001570 

0 

0.00471 

0 

3 

0.003141 

0 

0.00785 

0 

4 

0.004712 

0 

0.01099 

0 

5 

0.006283 

0 

0.01413 

0 

6 

0.007854 

0 

0.01727 

0 

7 

0.009424 

0 

0.02042 

0 

8 

0.01099 

0 

0.02356 

0 

9 

0.01256 

0 

0.02670 

0 

10 

0.01413 

0 

0.02984 

0 

11 

0.01463 

0.0026607 

0.02648 

6.801E-4 

12 

0.01259 

0.0123351 

0.02424 

1.337E-3 

13 

0.01183 

0 . 0197077 

0.02318 

1.923E-3 

14 

0.01141 

0.0266998 

0.02252 

2 .499E-3 

15 

0.01130 

0.0284328 
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After  defining  the  volume  and  area  of  the  control 
voliimes,  the  finite  difference  equations  ,as  shown  in 
equations  (A.  15)  to  (A.  17),  are  set  up  to  calculate 
temperature  of  each  layer.  The  time  step  (At)  is  chosen  as 
the  time  required  to  penetrate  through  the  first  layer. 

At  =  tpi  =  0.1578  Second 

The  heat  density  (HD)  of  each  layer  is: 
HDi  =  (piciVi,  i  +  P2C2V2,  i)Ah       for  i=l  to  15 
KAi  =  KiAi.  i  +  K2A2.  1  for  i=l  to  14 

For  the  first  layer  : 

Tl,  t  +  1  =   ^          Qin  +  (1  ^  ^)Ti,  t  +  — ^  

HDi  HDx  HDx 

For  the  layer  15 : 

At(fCAi4)  At(iCAi4), 

Tl5,  t  +  1  =    Tl4,  t  +  (1  )Tl5,  t 

HDis  HD15 

For  the  remaining  layers : 

Ti,  t  +  1  =   (Ti  -  1,  t  +  Ti  +  1,  t)  +  (1  ^  -yri.  t 

HDi  HDi 

for  i=2  to  14  (A.  17) 

After  setting  the  initial  temperatures  of  the  layers 
and  the  initial  current  time(tc)  to  zero,  the  finite 
difference  procedure  begins. 


(A. 12) 

(A. 13) 
(A. 14) 

(A. 15) 
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In  each  time  step:  tc  =  tc  +At 

For    tpk-i   <   tc    <    tpic,    applying    the    finite  difference 
equations    from   layers    1    to    k.       Due    the    theory   of  heat 
penetration  distance,    the   temperatures   of    the   rest  layers 
remain  the  same  as  the  temperature  of  the  last  time  frame. 
Ti,  t  +  1  =  Ti,  t         for     i  >  k  (A.  18) 

The  finial  temperatures  of  each  layer  can  be  determined 
after  the  current  time  reaches  the  final  heat  input  time 
(120  seconds) .  The  temperature  of  the  first  layer  represents 
the  temperature  of  the  position  of  the  heat  input. 


Table  A. 3 


Layer # 

Ave.  penetration 
distance  D 

Temperature ( °C ) 

1 

0.00 

240.73 

2 

0.01 

113 .88 

3 

0.02 

72  .83 

4 

0.03 

49.62 

5 

0.04 

34.46 

6 

0.05 

23  .98 

7 

0.06 

16.57 

8 

0.07 

11.31 

9 

0.08 

7.57 

10 

0.09 

4.93 

11 

0.10 

3.06 

12 

0.11 

1.92 

13 

0 . 12 

1.22 

14 

0.13 

0.76 

15 

0.14 

0.52 
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The  temperatures  of  the  rest  of  the  layers  represent 
the  temperatures  at  the  average  penetration  distances  (D)  of 
each  layer  as  shown  in  equation  (A.  19).  The  final 
temperature  and  average  penetration  distance  D  of  each  layer 
are  shown  in  Table  A . 3 . 

Di  =  0 

Dpi  -  1  +  Dpi  (A.  19) 

Di  =  —  — 

2 

The  temperature  distribution  throughout  the  beam  is 
obtained  by  expanding  the  temperatures  of  the  layers  to  two- 
dimensional  space.  First,  as  shown  in  figure  A.  4,  define 
the  origin  of  the  Cartesian  coordinate  system  at  the  left- 
end  bottom  of  the  beam  and  the  position  of  the  heat  input 
(0.3,0.2).  The  heat  penetration  distance  (Di)  to  a  point 
Pi(xi,yi)  in  material  1  can  be  calculated  from  equation 
(A.20) .  For  a  point  P2(x2,y2)  in  material  2,  the  penetration 
distance  from  this  point  to  the  position  of  the  heat  input 
can  be  calculated  in  equation  (A. 21)  from  the  theory  of  heat 
penetration  distance.  The  temperature  of  points  Pi  and  P2 
can  be  solved  as  shown  in  equation  (A. 22)  by  utilizing  the 
table  A. 3.  The  temperature  distribution  (T{x,y))  of  the 
beam  is  therefore  determined. 
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Di  =  VSa^03)^  +  (yr^0^  (A. 20) 


(A. 21) 


where  ri=hicsc(92) 

r2=V(x2  -  0.3)'  +  (ya  -  0.2)'  -  hiCscOz) 
e2=tan-Ml  (72-0.2) /(X2-0.  3)  I) 


For  any       >  D  >  Dk-i 

T(D)  =  T(Dk  -  i)  +  ( ^  ~       -  ^ )  (y(^)  _  y/Q^  _   N  (A.  22a) 

D;c  -  Dfc  -  1  ^  ^ 

for  and  D  >  0.14 


T(D)  =  0 


(A. 22b) 
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The  Thermal  Deflection  and  Stress  Distribution 

After  obtaining  the  temperature  distribution  (T(x,y)) 
of  the  beam,  the  neutral  axis  (shown  in  figure  A.  5)  of  this 
composite  beam  is  determined  by  utilizing  equation  (3.2) 
before  calculating  the  equivalent  moment. 

yo  +  O.l  yo+0.2 

j  EiydA  +  j  EiydA  =0  (A. 23) 

yo  yo  +  O.l 

From  which: 

Yo  =  -0.07473  m 


Yo 


T 


steel 


Aluminijm 


  Top  surface 

h,=0.1 

-*— Interface 


T 

Jh2=0 . 1 


•Bottom  surface 


Figure  A. 5  The  cross  section  of  the  beam 


The  total  flexural  rigidity  D  is  obtained  as  following: 

D  =  t[Ei{—  +  hi{—  +  yof]  +  k{E2(—  +  h2{hi  +  —  +  yo)'] 
12  2  12  2 

where  b:  width  of  the  beam  (A. 24) 

The  equivalent  moment  Mt(x)   is  obtained  as  follow: 

yo  +  hi  yo  +  hi  +  hi 

Mt(x)  =     j  EMiTix,  y)ydy  +     j  E2a2T(x.  y)ydy  (A.  25) 

yo  yo  +  hi 
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The  work  done  by  the  equivalent  moment  Mt(x)  is: 

L 

w  =  j  Mtv"dx  (A.  26) 

0 

By  assximing  the  deflection  function  v(x)    as  a  Fourier 


series ; 


V 


=   Ao+2^An  COS  (  )  +  2^  Bn  sin  (   )  ) 

L         „=i  L 


(A. 27) 


The  strain  energy  U  and  the  work  done  by  the  equivalent 
loading  (Mt  (x)  )  are: 


U  = 


Ein' 


^  (An^  +  Bn)n*  + 


n  =  l 
N  N 


X^-i                          "iiJ  ill 
>  {AnBn  —  —  (cos  (nji)  COS  (mji)  —  1) 


n  =  l  m  =  l 
m*n 


(A. 28) 


Ann^TI^J  Mt(x)  cos  {  )dx  -  2^  Bnn^Tt^J  Mt(x)  sm  (  )dx 

n  =  l  0  ^  n  =  l  0  ^ 

(A. 29) 

The  constraints  in  this  case: 
For  v=0  at  x=0 

Gi  =  Ao  +  ^  An  =  0 

n  =  l 

For  v'=0  at  x=0 


(A. 30a) 


G2  =  2"^^  -  0 
For  v"=0  at  x=L 

N 

G3  =  Ao  +  ^  An  cos  (2371)  =  0 

n  =  l 

For  V" ' =0  at  x=L 


(A. 30b) 


(A. 30c) 
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N 

Gi  =       Ann^  cos  (nTC)  =  0  (A.30d) 
ji=i 

For  v=0  at  x=Xs=0.7L 

nnxs^                    ,imxs^      „       ,     ^„  , 
Gs  =  Ao  +  2^  An  cos  {  )  +  +2^BnCOs(  )  =  0  (A.30e) 

The  Lagrangian (merit )   function  L  can  be  formulated  as: 

5 

L  =  U  -  w  +  ^  XiGi 

i=l 

By  equating  to  zero  the  partial  derivatives  of  the 
merit  function  with  respect  to  parameters  (1=1,2,3,4,5) 
Ao,  An  and  Bn  (n=l,  2,  .  .  .N)  ,  a  set  of  (2N+6)  linear  equations 
with  (2N+6)  unknowns  is  obtained.  For  N=8,10,14  the  unknown 
coefficients  of  the  deflection  functions  are  shown  in  Table 
A. 4.  The  results  from  the  proposed  method  for  different 
number  of  coefficients  N  are  shown  in  figure  A. 6. 


Table  A. 4 


N=8 

N=10 

N=14 

I 

sine 

cosine 

sine 

cosine 

sine 

cosine 

1 

-2  .92E-6 

1. 

92E-4 

-1.09E-4 

7. 

50E-3 

0.034467 

-1.10127 

2 

-3.37E-4 

9. 

50E-5 

-1.26E-2 

-2 

.56E-4 

1.92731 

0.12065 

3 

-1.55E-4 

-3 

.17E-4 

4.19E-4 

-1 

.39E-2 

-0.21661 

2.31035 

4 

2.46E-4 

-1 

.63E-4 

1.21E-2 

4. 

82E-5 

-2.24208 

-0.2765 

5 

1.32E-4 

1. 

42E-4 

-4.12E-4 

8. 

48E-3 

0.28756 

-1.85022 

6 

-6.04E-5 

7. 

62E-5 

-4.88E-3 

-2 

.78E-4 

1.32188 

0.24548 

7 

-3 .22E-5 

-1 

. 82E-5 

1.42E-4 

-2 

.26E-3 

-0.17696 

0.82149 

8 

1.66E-6 

-8 

.OlE-6 

8.08E-4 

5. 

31E-6 

-0.4426 

-0.10813 

9 

-1.29E-5 

2. 

05E-5 

0.05565 

-0.20461 

10 

-2.89E-5 

-7 

.15E-7 

0.07967 

0.02373 

11 

-0.00813 

0.02533 

12 

-0.00623 

-0.00212 

13 

3.77E-4 

-0.00107 

14 

9.88E-5 

3.42E-5 
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Figure  A. 6  The  results  from  the  proposed  method  for 
different  number  of  coefficients  N 


The  resultant  stress  is  the  superposition  of  the 
thermal  stress  (Ot)  ,  uniform  stress  (Cave)  and  bending 
stress (am).  In  order  to  satisfy  the  equilibrium  conditions, 
the  total  force  due  to  the  stress  distribution  in  x- 
direction  has  to  be  zero  and  the  total  moment  about  the 
neutral  axis  should  be  also  zero. 

For  material  1 : 


Ot  =  -aijE:iT(x,  y) 


El 


hiEi  +  hzEi 


^yo  +  hi  yo+hi  +  hj 

j    aiEiT{x.  y)dy  +      J  ttiSiTlx,  y)dy 

V.   yo  yo  +  hi 


(A. 31) 

^ 


(A. 32) 
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Cm  = 


Eiy 
D 


(  yo  +  M 


yo  +  /ii  yo  +  hi  +  h2 

j  ttiEiTU,  y)ydy  +      J  a2E2T(x,  y)Ydy 


where      yo  >  y  >  yo+hi 
For  material  2 : 
Gt  =  -(X2B2T(x,  y) 


yo+hi 


E2 


''yo  +  hi 


inEl   +  ^252 


y  0  +  hi  +  iij 


(A. 33) 
(A. 34) 


yo  +  hi 


J    ai£iT(x,  y)dy  +      J  aiBiT(x,  y)dy 

> 

(A. 35) 


B2y 


''yo  +  hi 


yo  +  hi+h2 


j  OCiEiT(x,  Y)YdY  +      I  CL2E^(x,  yfydy 


\^    yo  yo+hi 

where      yo+hi  >  y  >  yo+hi+h2 


(A. 36) 


The  resultant  stress  in  materials  1  and  2  of  the  beam 
can  be  calculated  as: 

O  =Ol+Oave+Om  (A. 37) 

The  maxim\am  temperature,  as  shown  in  table  A.  5  and 
figure  A. 7,  through  the  thickness  of  the  beam  occurs  at  the 
X  position  of  the  heat  input  (x=0.3m).  The  resultant  stress 
at  x=0.3  can  be  determined  by  using  equations  (A. 31)  to 
A.  37)  and  shown  in  figure  A.  8.  The  maximum  shear  stress  as 
shown  in  figure  A. 8  at  the  interface  of  the  two  materials  is 
60.9  Mpa. 
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Table  A. 5  The  temperature  distribution  through  the  thickness 


of  the  beam  at  x=0.3 


y-yo 

temperature 

y-yo 

temperature 

0.00 

240.73 

0.11 

2.59 

0.01 

113.88 

0.12 

2.12 

0.02 

72.83 

0.13 

1.75 

0 . 03 

49.62 

0.14 

1.46 

0.04 

34.46 

0.15 

1.19 

0.05 

23  .98 

0.16 

0.99 

0.06 

16.57 

0.17 

0.81 

0.07 

11.31 

0.18 

0.68 

0.08 

7.57 

0.19 

0.58 

0.09 

4.93 

0.20 

0.0 

0.10 

3.06 

Figure  A. 7  Temperature  distribution  in  the  beam 
at  x=0.3  t=120  sec 
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Slip  stress=60.9  Mpa 


ttttttttt'' 


<j>  <» 

to  to 

CM  tM 

00  r- 

O  T- 


« 

a 


Top 
surface 


interface 


Bottom 
surface 


Figure  A. 8  The  resultant  stress  and  the  shear  stress  between 

two  materials  at  x=0.3 
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